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The  Progress  Report 


This  report  is  the  forth  qucirtly  report  for  the  project  of  'Hocally  connected  adaptive 
Gabor  filter  for  real-time  motion  compensation,"  with  grant  number  N00014'94-l-0077, 
which  has  been  in  the  process  since  October  20th  of  1993  and  has  been  conducted 
under  the  supervision  of  the  principal  investigator,  Professor  Hua  Li  of  TexM  Tech 
University. 

As  at  the  end  of  the  first  year  of  this  three-year  project,  we  have  been  making  progress 
towards  the  goals  of  this  research  as  planned  in  the  project  proposal,  1(b)  on  page  21. 
In  peirticular,  we  have  been  accomplished  the  following  work  as  itemized  below: 

1.  We  have  conducted  extensive  reseMch  work  on  the  design  and  simulation  of 
electronic  analog  circuits  as  bcksic  building  blocks  for  the  VLSI  implementation. 
This  phase  of  the  work  is  a  little  bit  ahead  of  the  schedule  (about  1  month)  than 
we  originally  planned  in  the  proposal.  Tlus  hardware  design  phase  concurrent 
with  the  algorithm  analysis  and  verification  provided  coherent  work  <vnd  ensured 
the  quality  of  the  analog  VLSI  design.  The  work  in  VLSI  implementation  at  this 
stage  includes 

(a)  Design  one  of  the  most  essential  building  blocks,  a  video  frequency  opAmp. 
With  extensive  SPICE  simulations  using  the  device  model  provided  from 
MOSIS  actual  fabrication  run,  we  have  completed  the  design.  The 
characteristics  of  the  OpAmp  to  be  used  to  build  two-dimensional 
convolution  unit  is  given  in  a  research  memo  and  was  included  in  the  2nd 
quart  report  as  Appendix  II. 

(b)  Design  analog  multiplier,  which  is  to  be  used  for  building  two-dimensional 
analog  convolution  unit.  Ebctensive  SPICE  simulation  was  performed.  The 
preliminary  simulation  data  looks  very  positive  and  the  circuit  layout  design 
has  been  completed. 

(c)  Post-layout  SPICE  simulation  is  under  way  to  check  the  actual  VLSI 
implementation.  The  experimental  result  now  is  under  andysts  and  will  be 
documented  amd  reported  accordingly.  Some  results  have  been  documented 
and  given  in  Appendix  I  of  this  report. 

2.  In  order  to  benchmeu'k  the  VLSI  chips,  a  hardware  prototype  board  is  under 
design  and  construction.  This  board  will  be  used  to  compare  the  performance  of 
digital  approach  vs.  cuialog  approach,  analog  approach  bcised  on  the  standard 
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off-the  shelf  components  vs.  analog  customer-design  VLSI  approjrch.  The  board 
is  now  operating  and  it  is  tmder  calibration.  The  board  is  designed  with  a  486 
macliine  as  a  host.  The  prototype  demonstration  is  prepared  for  ONR. 

A  software  package  has  been  developed  to  implement  the  Gabor  motion  detection 
algorithm.  This  software  consists  of  utility  functions,  algorithm  modules,  and  test 
pattern  generators  for  the  experiments  and  verification  of  the  spatial  and  temporal 
selectivity.  A  software  program  manual  now  has  been  produced.  Fully  tested, 
documented  programs  are  provided  together  with  the  manual. 

This  software  package  consists  of  the  following  programs: 

1.  Programs  for  algorithm  implementation,  which  include: 

(a)  Real  Gabor  kernel  with  user-selectable  kernel  size,  spatial  frequency,  and 
orientation  frequency; 

(b)  Imaginary  Gabor  kernel  with  user-selectable  kernel  size,  spatial  frequency, 
and  orientation  frequency; 

(c)  Derivatives  (with  respect  to  x  and  y)  of  the  real  Gabor  kernel  with 
user-selectable  kernel  Mze,  spatial  frequency,  and  orientation  frequency; 

(d)  Derivatives  (with  respect  to  x  and  y)  of  the  imaginary  Gabor  kernel  with 
luer-selectable  kernel  size,  spatial  frequency,  and  orientation  frequency; 

(e)  Image  convolution  program  with  user-selectable  kernels; 

(f)  Least  square  estimation  algorithm  for  solving  an  over-determined  linear 
system,  which  is  needed  at  the  last  phase  for  optical  flow  computation. 

2.  Utility  programs  which  consist  of 

(a)  A  set  of  programs  to  creat  test  image  patterns  based  on  virtual  mlity 
technique.  These  programs  allow  user  to  define  the  observer's  position 
(camera  position),  the  fixation  direction  of  the  camera,  and  the  orientation 
of  the  projection  plane  (virtual  film.)  They  can  be  utilized  by  the  user  to 
define  objects  in  three-dimensional  world-coordinated  space. 

(b)  A  program  to  creat  animated  images.  This  program  will  produce  each 
individual  frame  of  images,  and  generate  an  animated  image  sequence. 
These  two  programs  will  be  particularly  needed  to  test  and  verify  the 
optical  flow  computation  based  on  the  created  known  motion  of  the  objects 
and  the  observer. 
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(c)  A  short  program  (in  MATLAB  format)  to  exeat  postscript  hies.  This  will 
produce  the  hardcopy  of  images,  and  computation  result. 

Based  on  the  theoretical  analysis,  the  optical  flow  was  computed  for  several  artiflcially 
generated  test  patterns.  These  test  patterns  are  designed  to  test  the  concept  of  spatial 
and  oritniation  selectivity.  The  patterns  were  generated  by  using  virtual  reality 
teclmique  based  on  three-dimensional  computer  graphics.  As  described  in  the  second 
quartly  report,  they  include  stationary  observer  while  object  is  moving,  and  stationary 
object  while  observer  is  moving,  as  well  as  both  observer  and  object  moving  in  a 
known  fashion.  In  order  to  further  test  the  computation,  we  have  modified  the 
resolution  and  color  depth  of  each  pixel  to  make  them  all  120-by-100  in  resolution  and 
S-bit  color  per  pixel.  A  set  of  floppy  disks  are  now  prepared  for  the  demonstration 
purpose.  This  set  of  floppy  disks  are  now  included  in  this  report  to  ONR.  Included  in 
this  set  of  disks  are  image  sequences  generated  under  known  conditions  as  benchmark 
for  the  future  testing  of  VLSI  implementation,  in  particular,  they  include  the  following: 

1.  An  observer  (camera)  is  stationary,  but  the  viewing  object  (a  sphere  with  radius 
equal  to  25  units)  is  moving  along  y-axis  5  units  per  frame. 

2.  The  observer  is  moving  along  y-axis  5  units  per  frame  while  the  viewing  object  is 
stationary. 

3.  Both  the  observer  and  the  viewing  object  are  moving  5  units  along  y-axis  and 
z-axis  respectively. 

As  briefly  summaried  here,  we  have  been  making  progress  as  planned.  In  the  next 
phase,  the  second  year  of  the  project,  we  will  produce  the  fabricated  analog  VLSI  chips 
for  testing  and  evaluation. 

The  End. 
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ANALOG  CONVOLUTION  UNIT  FOR 
REAL  TIME  IMAGE  PROCESSING 

Laszlo  Moldovan  and  Hua  Li 
Department  of  Computer  Science 
College  of  Engineering,  Texas  Tech  University 
Lubbock,  TX  79409 
E-mail:  xdhua@ttacsl.ttu.edu 


AhMraci —  The  dMign  of  an  analog  convolution  unit  for 
raal  tinM  Imago  procowlng  1*  proaontod  in  thlt  paper  with 
the  emphaals  on  the  deeign  of  the  batlc  building  block: 
CMOS  four  quadrant  multiplier.  Thlt  CMOS  multiplier 
operatet  with  and  bat  tingle  ended  voltage  input! 
making  it  very  eaty  to  uie  for  real  time  image  procetting. 
Thlt  bnaic  building  block  it  uied  to  build  a  5x8  array  M 
convolution  unit.  The  MAGIC  layout  of  tbit  circuit, 
ready  for  double  poly  CMOS  analog  fabrication  it  alto 
preaented. 

I.  INTRODUCTION 

Tlir  dmign  of  An  analog  convolution  unit,  for  real  time 
imago  procofwing  ia  given  in  tliia  memo.  Tlio  unit,  con- 
aiatj:  of  srvoral  <liifoi'ont.  kind  hnaic  liuit<ling  blorka  in¬ 
cluding  an  analog  multiplier.  An  analog  multiplier  can 
he  easily  huilt  in  bipolar  technology  and  they  have  been 
aueceiMfuliy  uae<|  for  yeara  |l).  The  problem  when  try¬ 
ing  to  tiiiihl  an  analog  multiplier  using  CMOS  devices  is 
tliat  tlic  output  current  of  the  sourec-rouple«l  tliffercn- 
tial  pair  (MO  and  M6  in  Pigiire  1)  drpen<ts  nonlinearly 
on  tlie  Idas  current  sinked  by  M7  and  the  input  voltage. 
Therefore,  the  linear  range  of  the  input  voll  ages  is  very 
limited  and  very  bard  to  eompensnte.  Given  the  nrni 
of  using  CMOS  tedinology  over  hipolar  tedinology.  we 
have  ilestgned  the  liasir  building  hloek  and  the  unit  liy 
using  artive  attenuator  tedinique.  Tliis  design  provides 
murli  better  linear  operating  range  of  the  CMOS  four 
qunilrant  multiplier. 

II.  PRINCIPLE  OP  OPERATION  OF 
ANALOG  MULTIPLIER 

Tlie  complete  circuit  of  the  multiplier  is  shown  in  Fig¬ 
ure  1.  The  core  of  this  circuit  consists  of  the  CMOS 
version  of  tlie  stanilard  Gilbert  cell  presented  in  [Ij.  As 
nl>ove  mentioned,  the  problem  of  this  cireiiit  is  the  re- 
ilurcrl  voltage  input  swing  range.  In  order  to  overcome 
this  problem.  I.iie  input  voltages  are  connected  to  the 
Gilbert,  cell  tliroiigb  artive  .itteiuiatars.  The  output  sig¬ 
nal  of  the  eireuit  is  a  dilferential  eiirrent.  wliirli  ran  he 
easily  convertwl  to  a  differential  voltage  with  two  load 
resistors  ronneetc«l  to  Vdd.  The  Gilbert  cell  and  the 
active  at.tenuators  are  ilew’rihed  as  follows. 

III.  GILBERT  CELL 

Devices  Ml.  M2.  M3.  Ml.  Mo.  MO  aiul  tlie  rurrent 
sink  M7  form  a  CMOS  vr-siim  of  the  Gilbert  roll.  All 


'^ransist.ors  operate  in  the  saturation  region  anti  tlie 
transconductancc  parameters  of  M1..M4  anti  M5,  M6 
are  matchctl  anti  equal  to  Ki  and  Jf  j  respectively.  There¬ 
fore,  the  iilcol  square-law  equation  can  he  iq>plietl.  Tlie 
output  currents  of  this  circuit  arc  I„i  anti  /oj.  given  by: 


loi  =  — (Lii  +  Ifo) 

(1) 

l<a  =  —[t,a  +  Lm) 

(2) 

Tims,  the  output  iliffcrential  eurrent.  Iad  —  f<a-  Ui  is 
given  by: 

/oa  =  ^/2^ -  .^1 

(3) 

where  Vosi  =  Vos4  =  V,'.  If  terms  anti 
are  mucli  smaller  than  1.  it  follows  that  /«/  tiepeniis 
Kncarly  on  V,  anil  is  given  liy: 

/«,  ~  (4) 

Also,  /ifi  anil  Ijt  ran  he  sliown  to  lie  in  tlie  following 
rdationsliip  witli  Feyj  s  Ven  =  V,: 

(5) 

Sulwtituting  (0)  into  (4).  wc  have  the  output  differ¬ 
ential  current: 

/„J  =  (C) 

wliidi  is  tlie  diaracteristir  of  an  analog  multiplier. 
But  tliis  equation  was  obtained  uniicr  tlie  assumption 
that  Vl  anti  ore  kept  very  small  wliidi  may  not  he 
satisfietl.  Tb  solve  this  problem,  these  two  voltages  were 
eolleetnl  at  the  outputs  of  active  attenuator^. 

rv.  ACTIVE  ATTENUATOR 

Tlie  active  attenuators  (one  for  the  VX  input  onotlier 
for  the  VY  input)  were  luiilt  with  transistors  M8.  M9 
anti  respectively  Mil.  M12.  In  ortler  to  couple  tlie  at¬ 
tenuators  with  the  two  inpiiU,  a  voltage  sltifl  was  ob¬ 
tained  with  tlie  source  followers  MIO.  MIS  anti  M13. 
M19.  A  voltage  (VX  or  VY)  npplicti  to  these  attenua¬ 
tors  plus  shifters  will  he  retlurnl  Iqr  a  factor  of  10.  as 
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resulted  from  PSPICE  simiilntions.  Tn  other  words,  the 
eireuit  will  rcjluee  the  slope  of  n  linearly  varying  voltage 
between  -IV  and  -|-1V  from  1  to  0.1. 

V.  MULTIPLIER  UNIT 

With  the  attenuators  attached,  t  he  transfer  character¬ 
istic  of  the  multiplier  unit  is  <lcscribe<l  by  the  equation: 

Id  =  hi  -  h2  ^  V2A'iA2  V.  (7) 

where  m  is  the  attenuation  factor  of  the  aeth’c  attenua¬ 
tors  (10.9).  So  far  the  output  was  a  «lifferential  current. 
If  we  want  to  convert  it  into  a  differential  voltage,  we 
nee<l  to  connect  two  matche<i  loa<l  resistors  between  the 
output.^  of  the  Gilbert  cell  and  V<l<l.  Thus  the  output 
voltage  is  given  by: 

Vo  =  /?/.  Id  ^  V.  Vy  (8) 

In  or<ler  to  be  able  t.o  ealibi'al.e  the  circuit  for  (lifferenl. 
fabrication  parainet.ors  (which  cannot,  he  known  pre¬ 
cisely  in  advance)  and  l.o  compensate  for  the  offset  volt¬ 
age.  the  atteinm.tors  iu'*'  biase<l  externally  on  the  gat.es 
of  M14  and  M15.  As  lomis.  DLl  amt  RL2  were  chosen 
to  be  IKrt.  The  sizes  of  the  l.raii'  l,.tors  are  given  in  the 
following  table; 


TABLE  I 

Device  sizes  for  multiplier  tinit 
size  W/L=/tin/pm 


li^ 

ma 

mm 

K'iHI 

iiisa 

Ena 

MIO 

mm\ 

ikm 

oa 

Mil 

mimi 

lixgl 

csoi 

mm 

ica 

Esm 

M13 

BQIl 

lEXTil 

msM 

Mil 

mu9sm\ 

lEa 

msM 

M15 

IDSai 

ica 

mm 
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VI.  EXPERIMENTAL  RESULTS 

First,  the  eireuit  w.v  test  ol  using  PSPICE.  The  transfer 
elinrneteristirs  for  a  -  IV  to  IV  de  sweep  on  VX  an«!  VY 
inputs  are  shfiwn  in  Figures  2  and  3.  The  PSPICE  file  of 
the  eireuit  is  cont.aine<l  in  APPENDIX  A,  A  Dotle  plot 
in  Figure  4  shows  that  the  circuit  has  a  -SdD  liandwidth 
of  24.4  MHz.  Following  the  preliminary  simulations,  the 
layout  of  this  circuit  Wiis  <lesigned  using  MAGIC  and 
then  its  circuit  was  extractcil  for  PSPICE  simulation 
to  verify  the  layout  design.  Tin-  PSPICE  siniulalions 
for  the  exte.-u-te*!  circuit,  gave  i.hr  t.r.ansfer  diaracteris- 
l.ies  its  shown  in  Figures  5  ami  0  and  the  Ou<le  plot, 
in  Figure  7.  The  MAGIC  layout  of  a  single  miiltiplier 
cell  c;ui  be  seen  in  Figure  8.  The  -3  <111  bantiwhllh  in 
this  <'i»se  is  19  MIlz.  Nonlim’iu  ity  tests  were  performetl 
on  tlie  an.tlog  multipli<'r  for  liotli  inputs.  First,  a  ±1V 
nuiip  signal  wiis  input  on  V'Y.  VX  w.-vs  ina«l''  succesively 
IV  iuid  -IV  ami  the  <iiitpnt  volt  itge  w.us  multiplied  with 


a  constant  factor  to  make  Vout/VY=l.  Then  VY  was 
subtractetl  from  Voiit  an<l  the  result  rcpresentctl  the 
absolute  nonlinearity  error.  FVom  this  value,  the  rela¬ 
tive  nonlinearity  error  was  calculated.  The  two  test  arc 
illustratetl  in  Figures  9  and  10  and  the  values  of  tlu; 
relat.ivc  errors  are  4.345%  an«i  respectively  3.61%.  Sim¬ 
ilarly,  the  same  signals  were  applietl  to  the  other  inputs 
and  the  some  operations  were  pcTformo<l.  For  this  cases 
we  have  the  situations  shown  in  Figures  11  and  12  an<l 
the  values  of  the  relative  nonlinearity  errors  are  4.438% 
and  respectively  0.12%.  In  order  to  find  the  total  har¬ 
monic  <fistor8ions  (THD)  of  this  circuit,  a  IMHz  sinu¬ 
soidal  signal  was  applied  to  one  of  the  inputs  and  the 
other  input  was  made  IV.  The  frequency  spectrum  of 
the  output  voltage  is  shown  in  Figure  13  and  the  THD 
due  to  the  second  and  third  harmonies  w.v:  ealeulaterl 
to  be  less  then  2%.  Using  the  stan<lar<l  40  pin  .analog 
frame  provi<le<l  by  MOSIS,  a  5x5  convolution  unit  an<l 
a  single  mtdtiplier  cell  were  combincxi  into  .a  single  mi¬ 
crochip.  The  5x5  convolution  unit  will  be  i.esteil  for  a 
T..a|>lncian  of  Gaitcsian  (LOG)  kernel.  The  kernel  values 
were  geiieratwl  using  the  LOG  formula: 


LOa(T.y) 


2tf*  - 1* 

2ff» 


if  I 


(9) 


The  values  for  the  kernel  were  calrulatetl  for  0=1  an<l 
are  sltawn  next: 


0.0519 

0.1231 

-0.1353 

-0.1231 

-0.0549 

•0.1231 

0.0000 

0.3033 

0.0000 

-0.1231 

-0.1353 

0.3033 

1.0000 

0.3033 

-0.1353 

-0.1231 

0.0000 

0.3033 

0.0000 

-0.1231 

•0.0549 

-0.1231 

-0.1353 

-0.1231 

-J.0549 

A  siinide  voltage  ilividrr  einniit  (Figure  14)  to  im- 
plctueiit  this  LOG  kernel  was  incorporated  in  the  mi- 
rrudiip  nn«l  ronueeted  to  the  Yij  inputs,  were  1  <  *  <  5 
«ud  1  <  j  <  5.  The  SPICE  file  of  the  voltage  <livider 
can  be  found  in  Appemlix  D.  All  the  Xij  inputs  were 
conneeteil  to  the  analog  input  pails.  The  load  resistors 
were  not  pul  in  endi  individual  cell.  Instead,  the  cor¬ 
responding  (/oL  and  Io2i.  1  <  *  <  25)  differential  cur¬ 
rent  outputs  were  conneeteil  together  to  provide  with 
the  sum  of  them.  Thus,  the  convolution  unit  performs 
also  the  siimmatinn  of  the  proilurts  of  the  25  multiplier 
cells  using  KCL,  The  result  of  the  coiivolutinn  can  be 
ai'cessoil  on  lol  and  Io2  pins  and  is  a  differenl.ial  current 
.  which  can  be  easily  roiivertcil  to  a  <lifter(>nti.al  voltage 
using  three  resistors.  Tliis  result  is.  therefore,  a  voll.age 
which  is  proportional  to  23*4 1  2Z*4i  XijY,,. 
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*TRANSFER  CHARARCTERISTIC  OF  THE  MULTIPLIER  FOR  A  VY  DC  SWEEP* 

Date/Time  run:  10/12/94  00:29:46  Temperature:  27. 
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BODE  PLOT  FOR  MULTIPLIER* 

Date/Time  run:  10/12/94  00:45:29  Temperature:  27. 


OHz  lOHz  lOOHz  l.OKHz  lOKHz  lOOKHz  l.OMHz  lOMHz  ICOMHz 

a  DB(I{Vo) /V(l) ) 

Frequency 


♦EXTRACTED  MULTIPLIER  TRANSFER  CHARACTERISTIC  FOR  A  VY  DC  SWEEP* 

Date/Time  run:  10/12/94  00:57:02  Temperature 


Figure 


.OliA  +  ‘ 


10/12/94 


ftaJ  \ 

W^si 

1  y^-s^  X  ^  tia 

k 

3 

■  ■■ 

— 

m 

n 

Ml  m  1 

:;;:i?iy;:::i5:x-i-- 

I 

■  ■  mm^ 

m 

1 

lP$S5ip;22si 

■  nB' 

1 

l^B 

Maa 

■  ■  K 

I 

■  m 

M 

•  ''  B 

m  bH 

B  m 

Date/Time  run:  10/15/94  11:96:38 


•NONLINEARITY  ERRORS  OF  ANALOG  HULTIPLIER* 
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A V<3,5) *672.45-  V(13) 


54-  V(10) 


NONLINEARITY  ERRORS  OF  ANALOG  MULTIPLIER* 


APPENDIX  A 


ANALOG  MULTIPLIER  WITH  ASYMMETRIC  INPUTS 


.OPTION  NOECHO  NOMO.O 

♦DEFINITION  OF  MODELS 

♦N43B  SPICE  LEVEL  2  PARAMETERS 

.MODEL  N  NMOS  LEVEL-2  PHI-0.600000  TOX-4.2600E-08  XJ-0.200000U  TPG-1 
+  VTO-0.8109  DELTA-6.7500E+00  LD-9.2070E-08  KP-4.8109E-05 
+  UO-593.5  UEXP-1.4990E-01  UCRIT-8.1140E+04  RSH-2.4540E+01 
+  GAMMA-0.4512  NSUB-4 . 0300E-t-15  NFS-1. 98E-(-ll  VMAX-6.2580B-t-04 
+  LAMBDA-2. 62 lOE-02  CGDO-1 . 1195E-10  CGSO-1.1195E-10 
+  CGBO-4.7024E-10  CJ-1 . 0922E-04  MJ-0.8608  CJSH-2 .4755E-10 
+  MJSW-0. 035834  PB-0. 800000 

*  Weft  -  Wdrawn  -  Delta_W 

*  The  suggested  Delta_W  is  -7.8340E-07 

.MODEL  P  PMOS  LEVBL-2  PHI-0.600000  TOX-4 .2600E-08  XJ-0.200000U  TPG— 1 
+  VTO— 1.0600  DELTA-9. 0900E+00  LD-1 .7890E-07  KP-2 .2267E-05 
+  UO-274.7  UEXP-3.4430E-01  UCRIT-6 . 9200E+04  RSH-5 . 9420E+01 
+  GAMMA-0.3746  NSUB-2 . 7770E+15  NFS-3. 23E+11  VMAX-9. 9990E+05 
+  LAMBDA-5. 3 9 90E-02  CGDO-2 . 1752E-10  CGSO-2 . 1752E-10 
+  CGBO-3.9953E-10  CJ-3 . 1585E-04  MJ-0.5914  CaSW-3.2974E-10 
+  MJSW-0. 315516  PB-0. 700000 

*  Weff  -  Wdrawn  -  Delta_W 

*  The  suggested  Delta_w  is  “3.4000E-07 

*ln,**lni1i**it**»********»**********************iiii************it1t*-kit*1i******** 


♦GILBERT  CELL 


Ml  3  2  1  4  N  M-80U  L-4U 
M2  5  17  1  4  N  W-80U  L-4U 
M3  3  17  6  4  N  W-80U  L-4U 
M4  5  2  6  4  N  W-80U  L-4U 
M5  1  18  7  4  N  W-80U  L-4U 
M6  6  8  7  4  N  W-80U  L-4U 
M7  7  19  4  4  N  W-80U  L-4U 


♦ATTENUATOR  1 


M8  11  10  9  9  P  W-2U  L-3U 
M9  4  10  11  9  P  W=2U  L-8U 
MIO  9  11  2  4  N  W-2U  L-3U 
M14  2  16  4  4  N  W-16U  L-2U 


♦ATTENUATOR  2 


Mil  12  13  9  9  P  W-2U  L-3U 
M12  4  13  12  9  P  W-2U  L-8U 
Ml 3  9  12  8  4  N  W-2U  L-3U 

M15  8  16  4  4  N  W-16U  L-2U 

RLl  3  9  IK 

RL2  5  9  IK 

RL  3  5  lOMEG 


Vdd  9  0  5V 
Vss  4  0  -5V 
Vbias  16  0  {VBl 


V17  17  0  OV 
VIS  13  0  OV 


V19  19  0  OV 


*VX  10  0  {X) 

VX  10  0  SIN<0  1  IMEG) 

*VX  10  0  AC  IV 

*VX  10  0  PWL(0US  {-XJ  2US  {X}  4US  {-X>  6US  {X)  8US  {-X)  lOUS  {X}) 

*VX  10  0  PWL(ONS  OV  200NS  OV  200.  COINS  IV  600MS  IV  600.001NS  OV  lUS  OV) 
VY  13  0  {Y} 

*VY  13  0  SIN(0  1  IMEG) 

*VY  13  0  AC  IV 

*VY  13  0  PWL(ONS  {-Y)  2US  {Y}  4US  {-Y}  6US  {Y)  8US  <-Y)  lOUS  {Y}) 

.PARAM  X=1V 
.PARAM  Y»1V 
.PARAM  VB— 3.915 

♦.STEP  PARAM  Y  -1  1  .1 
*.DC  PARAM  Y  -1  1  .1 
.TRAN/OP  .OIPS  lUS 
*.AC  DEC  100  1  1E8 
♦.WATCH  DC 
.OP 

♦.PROBE  V(14)  V(18) 

.END 


r 
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♦voltage  dividers  for  kernel  coefficients* 

■OPTION  NOECHO  NOMOD 

♦DEFINITION  OF  MODELS 

♦N43B  SPICE  LEVEL  2  PARAMETERS 

■MODEL  N  NMOS  LEVEL-2  PHI-0.600000  TOX-4 .2600E-08  XJ-0.200000U  TPG-1 
+  VTO-0.8109  DELTA-6. 7 500E+00  IjD=9 . 2070E-08  KP-4 . 8109E-05 
+  UO-593.5  UEXP-1.4990E-01  UCRIT-8 . 1140E+04  RSH-2 . 4540E+01 
+  GAMMA-0.4512  NSUB-4 , 0300E+15  NFS-1. 98E+11  VMAX-6.2580E+04 
+  LAMBDA-2. 62 lOE-02  CGDO-1 . 1195E-10  CGSO-1 . 1195E-10 
+  CGBO-4.7024E-10  CJ-1 . 0922E-04  MJ-0.8608  CJSW-2 .4755E-10 
+  MJSW-0. 035834  PB-0. 800000 

*  Weff  =  Wdrawn  -  Delta_w 

*  The  suggested  Delta_W  is  -7.8340E-07 

■MODEL  P  PMOS  LEVEL-2  PHI-0.600000  TOX-4 .2600E-08  XJ-0.200000U  TPG— 1 
+  VTO— 1.0600  DELTA-9. 0900E+00  LD-1 ■7890E-07  KP-2 . ^267E-05 
+  UO-274.7  UEXP-3.4430E-01  UCRIT-6 . 9200E+04  RSH-5 . 9420E+01 
+  GAMMA-0.3746  NSUB-2 . 7770E+15  NFS-3. 23E+11  VMAX-9 . 9990E+05 
+  LAMBDA-5. 3 990E-02  CGDO-2 . 1752E-10  CGSO-2 . 1752E--10 
+  CGBO-3.9953E-10  CJ-3 . 1585B-04  MJ-0.5914  CJSW-3 .2974E-10 
+  MJSW-0. 315516  PB-0. 700000 

*  Weff  -  Wdrawn  -  Delta_W 

*  The  suggested  Delta_W  is  -3.4000E-07 

*ir***lr************ir*****ir***1i1ricir*****1r*1tlr*1tik****it*1r**********it*it1tit1r 


♦KERNEL 


♦for  IV 

M16  14  22  4  4  N  W-lOU  L-2U 
M17  15  15  14  4  N  W-2U  I.-6U 
M18  9  9  15  4  N  W-3U  L-2U 

♦for  0.3033V 

Ml 9  20  22  4  4  N  W-8U  L-2U 
M20  21  21  20  4  N  W-6U  L-2U 
M21  9  9  21  4  N  W-2U  L-2U 

♦for  -0.055V 

M22  23  22  4  4  N  W-8U  L-2U 
M23  24  24  2J  4  N  W-7U  L-2U 
M24  9  9  24  4  N  W-2U  L-3U 

♦for  -0.123V 

M25  25  22  4  4  N  W-7U  L-2U 
M26  26  26  25  4  N  W-2U  L-6U 
M27  9  9  26  4  N  W-2U  L-4U 

♦for  -0.1353V 

M28  27  22  4  4  N  W-7U  L-2U 
M29  28  28  27  4  N  W-2U  L-4U 
M30  9  9  28  4  N  W-2U  L-4U 


**************************************^k*************** 


VDD  9  0  5V 
VSS  4  0  --5V 


Vbias  22  0  {VB) 


.PARAM  VB— 3.5V 

•■STEP  PARAM  VB  -3.6  -3.4  .01 

*.DC  PARAM  VB  -4  -3  .1 

.TRAN  IPS  2NS 

*. PROBE 

.OP 

.END 
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PREFACE 


This  software  manual  <lcscTihcs  the  algorithms  that  have  been  «lcvclopc<l  for 
the  <lemonstration  project  on  motion  detection  and  motion  compensation  using 
biologically  inspirc<l  Gabor  transforms.  The  objectives  of  the  project  arc  to  (1) 
<lcvelop  an  algorithm  suital)le  for  real  time  analog  VLSI  (Very  Large  Seale 
Integi  atc<l  Electronics  Circuits)  implementation  and  (2)  fabricate  the  design 
through  MOSIS.  The  algorithms  <levelopc<l  and  teste<l  in  this  package  will  be 
used  as  bench  marks  to  objectivdy  evaluate  the  performance  of  the  VLSI  chips, 
and  the  prototype  boar<i.  This  package  is  a  part  of  the  result  of  the  rescarcli 
project  support.c<l  by  ONR> 

This  package  inclu<lc8  three  categories  of  programs:  (1)  The  programs  that  are 
usc<l  to  generate  test  pattern  images.  (2)  The  programs  that  are  developed  to 
compute  optical  flows  by  using  Gabor  functions,  or  Gabor  Transforms.  And  (3) 
the  programs  that  arc  use*!  to  generate  sequence  of  images  by  using  virtual  re¬ 
ality  and  tlirtxsclimcnsional  computer  ^aphics  tccluiiqucs.  During  the  process 
of  preparing  this  package,  we  have  extensively  teste*!  eacli  program.  Most  of 
the  programs  in  this  package  were  written  in  ANSI  C  to  ensure  the  portability 
across  <lifFcrcnt  hardymre  platforms.  Fbr  the  programs  to  generate  3D  virtual 
reality  studio  environment,  we  have  used  the  program  language  <levclopc<l  by 
Watkins  which  is  a  C-type  language.  Althougli  most  of  tiic  programs  were 
<lcvclope*l  for  SUN  SPARC  Station,  they  can  also  be  compile*!  for  486  personal 
computer  with  some  minor  mo*liflcation8.  The  use  of  these  programs  for  ONR 
is  free,  but  we,  the  authors  of  these  programs  an*l  manual,  make  no  warranty  of 
any  kiml,  expresse*!  or  implicsl,  with  regard  to  the  programs  or  the  documen¬ 
tation.  Therc^forc,  the  authors  shall  not  be  liable  in  any  event  for  incielcntal 
or  consequential  damages  in  connection  with,  or  arising  out  of,  the  furnishing, 
performance,  or  use  of  these  programs.  All  branel  names,  trailemarks,  arc  the 
property  of  their  respesetive  holders. 

Several  researcJi  assistants  at  Texas  Tech  University  arc  responsible  for  *1*> 
vcloping,  implementing,  an*l  testing  of  the  programs.  The  present  version  of 

C.D.  Watkinii,  S.D.  Coy,  nnd  M.  Finlay,  PhntortaliDm  nn4  Hay  IVneing  in  C,  M&T 
nook*.  New  York,  NY  lOOll,  1992, 


iii 


IV 


the  manual  was  proof  read  and  tested  at  la's  Latkoratory,  Computer  Science 
Department,  College  of  Engineering,  Texas  Tbdi  University  the  prindpal 
investigator  of  the  project.  Inquiry  an<i  questions  regarding  the  source  code, 
the  algorithms  sliouUl  he  <lirccted  to  the  principal  investigator. 


Principal  Investigator:  Hua  Harry  Li,  Ph.D. 

Computer  Science  Department 

Collgc  of  Engineering 

Luhhock,  Texas  79409 

Phone:  (806)  742-3513 

Fox:  (806)  742-3519 

E-mail:  huaOcs.ttu.c<lu 
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INTRODUCTION 


ABSTRACT 

A  brief  summary  of  the  programs  and  tlieir  functionalities  included  in  this  manual 
are  (;iveu  in  this  chapter.  Namely,  three  major  categories  of  programs  are  descrilwd. 
These  categories  are;  (1)  the  programs  for  generating  test  pattern  images  as  liencli 
marks,  (2)  the  programs  for  computing  optical  flow  basol  on  Galtor  functions,  or  Ga- 
lior  transforms,  and  (il)  the  programs  base*!  on  virtual  reality  teclinique  for  generating 
3D  camera  mo<lcls  sumI  creating  user-defined  and  user-controlled  image  patterns. 


1  BRIEF  OVERVIEW 

Tliis  package  includes  three  categories  of  programs: 

1.  The  programs  that  are  used  to  generate  test  pattern  images.  Tlic  images 
are  ciiaracterized  with  uscr-<lefincd  orientation  and  spatial  frequence. 

2.  The  programs  that  are  <lcvclopc<l  to  compute  optical  flows  by  using  Ga¬ 
bor  functions,  or  Gabor  IVansforms.  In  this  category,  there  are  utility 
programs  for  solving  a  set  of  linear  algebraic  system  by  using  least  square 
estimation  (LSE)  tccliniques,  and  programs  for  user  selectable,  user  mo«l- 
iflable  Gabor  kernels. 

3.  The  programs  that  are  used  to  generate  sequence  of  images  by  using  virtual 
reality  and  thrce-<limcn8ional  computer  graphics  techniques.  Users  arc 
given  the  control  of  selecting  the  optical  chttracteristics  of  a  virtual  camera 


0 


6 


Chapter  1 


to  generate  wcll-controllcxl  eaincra  motion,  tw  ego  motion.  This  then  is 
uscil  to  precisely  cliorncterizc  the  amount  of  self  motion  in  comparison  to 
object  motion. 


In  this  <locumcntation.  we  use  the  term  “two-dimensional  motion"  to  deserihe 
the  image  pattern  generated  <Hrectly  by  shifting  image  plane,  which  does  not 
involve  the  selection  of  camera  mo<Iel,  the  calculation  of  perspective  projection. 
While  the  term  “thrce-<limcnsional  motion"  referres  to  as  the  image  patterns 
generated  by  using  virtual  reality  technique.  Namely,  the  teclmiquc  <lefines 
camera  characteristics  an<l  creats  realistic  looking  perspective  projection  im¬ 
ages.  The  two-dimensional  motion  patterns  are  majorly  used  to  calibrate  the 
orientation  selectivity  and  spatial  frequency  selectivity  of  the  Gabor  kernels. 


2  PREPARING  PROGRAM  FOR 
EXECUTION 

The  compiled  and  linkc<l  executables  of  the  programs  are  provi<lc«l  on  a  floppy 
disk  (MSDOS  readable  format)  with  this  manual.  These  executables  are  pre¬ 
pares!  under  the  following  c.on<litions; 


1.  The  operating  system:  SunOS  4.1.2. 

2.  The  cc  compiler. 

3.  The  hardware  platform;  SUN  SPAI1.C  1  workstation. 


If  you  have  the  complieel  an<l  linked  executables,  in  order  to  run  these  programs 
all  you  have  to  do  is  to  type  the  program  name  and  hit  return.  If  you  only 
have  the  source  c.o<lc  and  yuu  need  to  compile  and  link  them,  then  follow  the 
instructions  given  below: 

cc  -o  output-filename  inputfile.e  -Im 


whicli  will  creat  the  executable. 
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GENERATING  2D  TEST  PATTERN 

IMAGES 


ABSTRACT 

Two  difTeMiit  types  of  test  pattern  images,  3D  test  patterns,  and  3D  test  patterns, 
have  lieen  utiliietl  for  testing  the  motion  detection  algoritlun  as  explaine<l  in  Section 
1.2.  The  2D  test  patterns  are  grey  level  sine,  cosine,  or  strip  hw  images.  These 
test  patterns  are  use<l  to  csdihrate,  test  stnd  demonstrate  the  orientation  and  spatial 
frequency  selectivity  of  Gabor  kernels. 


1  CREATING  COSINE  WAVE  IMAGE 
PATTERNS 

DECRIPTION;  crcat  a  2-D  cosine  wave  image. 

USAGE:  Co8lmg(int  row,  int  column,  float  omega,  float  orient,  float  spcc<l, 
float  time) 

int  row:  row  of  image  in  pixels. 

int  column  :  column  of  image  in  pixels. 

float  omega:  spati^  frequency  of  image  in  cycles  per  pixel. 

float  orient;  orientation  of  image  in  degree. 

float  spectl:  the  variable  to  define  tl»c  phase  sliift. 

float  time;  the  vnrinl)lc  for  defining  the  phase  sliift. 

REMARKS: 
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Chapter  2 


Coaling  generates  a  2*D  cosine  wave  imogc  witli  user  sdcetnl  resolu¬ 
tion  (row-lty-eolumn)  and  user  dcfincti  spatial  frequency  and  orienta¬ 
tion  of  the  cosine  pattern. 


MATHEMATICAL  FORMULATION: 

y)  =  co«(2  ♦  p»  ♦  otnega  *  (x  *  cx>8(Vu:ta)  +  y  ♦  «sn(t/tcta)  -  speed  *  time)). 

where  variable  omega  defines  a  spatial  frequency  in  terms  of  cycles  per  pixel, 
variable  theta  ileRnes  on  orientation  of  the  test,  pattern.  The  variolile  spceil  de¬ 
fines  the  phase  shift  of  the  image  pattern  wiiirh  can  be  interpretcii  as  a  starting 
time  instance  when  the  test  pattern  image  was  creatcil.  Dy  clioosing  different 
pliasc  shift,  ii.^er  con  creat.  a  sequence  of  images  arranged  in  tlic  order  of  small¬ 
est  phase  shift  to  the  biggest  phase  shift  to  represent  a  sequence  of  moving 
cosine  patterns.  The  variolile  time,  tc^ether  with  the  varialilc  speed  creats  the 
desired  phase  sliift. 


EXAMPLE: 


iincluds<stdllb> 

*includs<stdlo  .li> 

Int  mainO 

{ 

char  *fiI«nan««"outpat.dat"; 
int  rov*31,  coluiiin=31; 

float  spatial.f raqusneys .  1 ,  oriastatloBi^SO : 
float  spsad  »  i,  time  »  0.4; 

CosImg(rov,  column,  spatial.frtqusncy,  oriantatloa,  spasd,  tlmt); 
rsturn  0; 

> 


By  chenging  the  different  phase  shift  value,  the  sine  wave  function  can  be 
createil  using  the  same  program. 


2  CREATING  BAR  STRIPE  IMAGES 


DEscnii*  rioN 


USA  (i  10:  I]:i.rlini'.!  ini  ri  iw.  ini  •■ilninn.  I  Ion  I.  onnT,;i.  Ilo.i.i.  oriim..  rlmr’  lilcna 


REM  A  IJKS 
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M ATH M ATIC A  L  FORM  IJ  L ATIO N : 

f  255  if  <-iis27r(.fu;i  +  >  0.0 

II  ol.liorwi.M* 


Fimir*!  '2.  A  l>;ir  Ht.rii>  witli  oim^git  cqiiiil  fo  (1.05  BHii  thntJi 

(V]U»1  t.<>  (I  <(i!(!;r<xv 


EXAMPLE; 


#include<stdlib> 

#include<stdio .h> 

int  mainO 

{ 

char  ♦f  ilenaine="output  .dat" ; 
int  row=.'51 ,  column=31 ; 

float  spatial  .f  requency=  .  1  ,  or  i  eiitat  ion=.30 ; 

BarTingfrow,  column,  spatial  .frequency ,  orientation,  filename); 
return  C ; 

) 


_ 3 

MOTION  DETECTION  BASED  ON 

GABOR  FUNCTIONS 


ABSTRACT 

Tiiis  cliapter  th«  program*  that  implement  motion  iletection  algoritlim  l)a«eil 
on  Galior  functions.  GaUor  functions  are  very  special  kind  of  functions  wliicli  mimics 
orientation  selective,  spatial  frequency  selective  characteristics  of  human  early  vision 
system.  Included  here  are  the  functions  for  creating  Gahor  kenels,  the  utility  func¬ 
tions  for  computing  convolution  of  real  anil  complex  Gabor  kernels,  the  fimctions  for 
solving  a  set  of  linear  algebraic  system  baseii  on  least  square  estimation  (LBS),  and 
tiie  function  for  computing  motion  vectors,  as  well  as  the  utility  functions  for  ploting 
optical  flow  pattern. 


1  CREATING  GABOR  KERNELS 

DECRIPTION:  Great  2-D  Gabor  kernels  ami  the  rlcrivatives  (with  respect 
to  X  and  with  respect  to  y)  of  Gabor  kernels. 

USAGE: 

void  CrcatcGabor(oniega,  theta,  type,  size,  tempt,  tempx,  tempy) 


float  omega:  modulateil  spat.ial  frequency  (uo)  (cycles  per  pixel). 

float  theta:  orientation  of  the  kernel  (degree). 

int  type:  kernel  type:  1  for  real  kernel,  0  for  imaginary  kernel, 

int  size:  the  rlcsired  kernel  size. 

float**  tempt:  pointer  of  a  block  to  Gabor  kernel. 
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float**  tcmpx:  pointer  of  a  block  to  the  derivative  (with  respect  to 

x)  of  the  Gabor  kernel. 

float**  tompy:  pointer  of  a  block  to  tlic  derivative  (with  respect  to 

y)  of  the  Gal>or  kernel. 

REMARKS: 


CrcatcGaborO  ercats  Gabor  kernel  and  its  derivative  kernels  (the 
derivative  with  respect  to  x  or  y.)  The  moduh'ited  spatial  frequency, 
orientation,  and  size  of  the  kernels  arc  determined  by  tlie  user  8clccte<l 
omega,  theta,  and  size.  Tlic  result  of  the  created  kernel  is  stored  in 
three  blocks  dynamically  allocated  in  the  memory. 

MATHEMATICAL  FORMULATION: 

Gr(x,v)  =  c*  ^Ws2x(w,a  + 

y)~c^~  *  •'**  ^«n2ir(w,a!  +  UyP). 
the  derivatives  of  Gatior  real  and  imaginary  parts  arc 

=  —2tcc  7^(^  cos  27r{wsX  +  UyV)  +  w,  sin  +  Uyp)), 

=  -2»ro  ^(^co8  2w(w,3  +  Wyy)  +  u;ySin2ff(a;,®  +  Wyv))T 

OG-  -sldUi  X 

=  —2kc  ( - j  sin  2w(w,3;  +  w^y)  +  w,  cos 2n'(b;aX  +  UyV)), 

vX  O* 

=  — 27rc“*  (~^  sin  2»r(u),*  +  w^y)  +  Wy  cos 2?r(b;aZ  +  UyU))- 

where  w,  and  uiy  arc  mo<lulatcd  spatial  frequency  components  in  sptial  fre¬ 
quency  domain, 

EXAMPLE: 


#includ«<stdlib> 
#lnclude<8tdio .h> 
#iacludG<itiath  .h> 
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int  mainO 

•C 

static  float  ••k«rn_t,  ••kam_y; 

float  spatial_fraq  =  .1.  oxlantation  =  0; 
int  filter_type  =  0; 
int  TOW,  col,  kernsiza,  indax; 

FILE  aoutt,  aoutx,  ♦outy; 

/*  calculata  kemal  size  and  indax  */ 
kernsize  =  (int)cail(l .0/(spatlal_fraq)) ; 

if (f abs ( (f loat)keTn8iza/2 .0-cail( (float )karnslza/2 . 0) ) <=0 . 01 ) 
karnslze-i"*-; 

index=(int} (((float)kernsiza-1.0)/2.0) ; 

/*  allocate  memory  for  kernels  */ 
kern.t  s(fioat  **)  mallocCkemsizeesizeof (float*)) ; 
kern.x  =(float  **)  malloc(kemsize*slzaof (float*)) ; 
kern_y  ‘^(float  **)  malloc(kemsize*8izeof (float*)) ; 
for  (1*0:  i<  kernsize :  i++  )  { 
kern_t[i]  =  (float  *)malloc(kem8ize*sizaof  (f loat)*l)  ; 
kern.xCl]  =  (float  •)malloc(kam8ize*Blzaof (float)*l) ; 
kern_y[l]  =  (float  *)malloc(kem8lze*Bizaof (float)-*-!) ; 

> 

/*  generate  Gabor  type  kernels  */ 

CreateGabor(Bpatial_freq,  orientation,  fllter_type,  index,  kem_t,  kem_x,  kern.y) ; 

if((outt  =  fopen("gabor_t.dat",  "w"))  ==  NOLL)  exit(l); 
if((outx  =  fopen("gabor_x.dat",  "w"))  ==  NULL)  exlt(l); 
if((outy  =  fopen("gabor_y .dat",  "w"))  ==  NULL)  exlt(l); 
for(i=0|  i<  kernsize:  i-*-*-)  i 
for  (j=0:  j<  kernsize:  j-f-*-)  { 

fprintf(outx,  "%f\t",  kem_t[i]  [j]) : 
fprintf  (outx,  "j(f\t",  kem_t[i]  [j]) : 
fprintf  (outx,  '"/,f\t",  kem.t  ti]  [j] ) : 

> 

fprintf (outx,  "Nn"): 
fprintf (outx,  "\n"); 
fprintf (outx ,  "\n"): 

} 


f close (outt) : 
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f close (outx) ; 
f close (outy) ; 
return  0; 

} 


2  COMPUTING  CONVOLUTION  OF 
GABOR  KERNELS 

2.1  Image  Buffer  Initialization 

DECRlPTIONs  initialize  a  buffer  for  convolution. 

USAGE:  voi<l  Duffcrlnit(in<lcx,  kemsizc,  column,  buffer,  fp) 


int  imlex:  index  of  Gabor  type  kernel. 

int  kernsize;  Gabor  kernel  size 

int  column:  column  of  the  given  image. 

unsigned  char**  buffer:  pointer  of  the  block  storing  image. 

FILE*  fp:  pointer  of  an  image  Blc. 


REMABICS: 


Bufferlnit()  initializes  a  buffer  for  convolution.  The  use  of  a  buffer 
accomplishes  the  convolution  of  a  given  image  one  row  at  a  time.  The 
size  of  the  buffer  can  be  dcterminecl  by  the  product  of  the  number 
of  rows  of  a  given  kernel  and  the  number  of  columns  of  the  given 
image.  The  top  half  space  of  the  buffer  is  initialized  to  be  filled  with 
the  first  row  of  the  given  image  for  the  consideration  of  convolution 
with  image  boundary  and  another  half  of  the  buffer  was  filled  with 
the  image  <lata.  The  image  <lata  sliouhl  be  8-bit  per  pixel  in  binary 
format.  The  buffer  is  a  static.  unsigne<l  char  double  pointer  type. 


EXAMPLE: 


#lnclude<8tdlib.h> 
#lnclude<stdlo  .li> 
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#l&clud«<malloc  .h> 
lat  mainO 

static  iiuslgnad  char  **lmagabufl«r; 

Int  image _row=50,  image _col=EO: 

int  keriiel_i&dex=7.  kemel.size  =  IS; 

int  i,  j; 

FILE  ein; 

/•  allocate  memory  for  image  buffer  */ 

imagebuffer  =(unsigiied  char  **)  malloc(image_rovesizeof (unsigned  char*)); 
for  (i=0;  i<  image.col;  i++  )  { 

imagebuffer [i]  =  (unsigned  char  *)malloc(image_col*sizeof (unsigned  char)); 

> 

/*  open  image  file  in  READ.ONLY  and  BINARY  mode  */ 
if ((in  a  fopenC'image.dat",  "rb"))  ==  NULL)  exit(l); 

/*  initialize  image  buffer  «/ 

Bufferlnlt(kemel_index,  kemel.size,  image.col,  imagebuffer,  in); 

ford  =  0;  i  <  kemel.size;  l-«^  ) 
for  (j  =0;  j  <  Image.col;  j  ++  ) 

printf ("Image  Data  tXdlCXd]  ■*  Xd\n",  i,  j,  imagebuffer [1]  [j]) ; 

return  0; 

} 


2.2  Image  Buffer  Shifter 

DECRIPTION:  up<latc  the  content  of  the  image  buffer  sliifting  up  one 
row.  This  shift  will  result  row  1  of  the  buffer  being  iliscardcd,  row  2  movctl  to 
row  1,  row  3  movc<l  to  row  2,  etc.,  an<l  a  new  image  row  movefl  to  the  last  row 
of  the  buffer.  The  shift  is  pcrformc«l  every  time  tiic  convolution  of  one  row 
finished. 


USAGE: 

void  BufferShifter(kernjndcx,  kem.sizc,  rowJmg,  colJmg,  rowJoop,  buffer,  fp) 
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int  kcmjmlcx:  index  of  Galior  kernel, 
int  kcm^zc:  G&ltor  kernel  size, 
int  rowjmg:  row  of  image, 
int  columg:  column  of  image. 

int  rowjoop;  actual  row  number  in  convolution  double  loop, 
unsigned  char**  buffer:  pointer  of  a  block  storing  image. 
FILE*  fp:  pointer  of  the  image  file. 

REMARKS; 


BuflFcrShiftcr()  shifts  <lata  of  each  row  in  image  buffer  up  to  next  row 
an<i  up<latcs  data  of  last  row  in  image  buffer  after  each  convolution 
operation.  The  rowJoop  is  the  actual  row  number  in  convolution 
<loublc  loop  in  onlcr  to  consider  the  lower  boundary  in  the  convolution 
operation.  The  image  format  should  be  8-bit  per  pixel  in  binary,  and 
the  image  buffer  could  be  static.  un8ignc<l  char  double  pointer  type. 

EXAMPLE; 


*inclu<le<8tdlib  .h> 

#includ«<stdio  .}i> 

#'  lude<malloc.h> 

int  mainO 

i 

static  unsigned  char  eeimagabuffsT; 
int  imags.xovsSO,  imaga_col=60 ; 
int  kexn«l_indax=7,  kamal.siza  =  15; 
int  i,  j; 

FILE  *in; 

/*  allo'  msBoxy  lor  inage  buffer  •/ 

:  igeb'.'  -(luisigned  char  «*)  oalloc(image_rov*slzeof (unsigned  chare)); 
for  (i=0;  i<  image.col;  i++  )  { 

image buff ex [i]  =  (unsigned  char  *)nalloc(image_cole8lzeof (unsigned  char)); 

> 

/*  open  itt  ■  file  in  READ.ONLY  and  BINARY  mode  »/ 
if  ((in  =  ;  .^^-uC'image.dat",  "rb"))  ==  NULL)  exlt(l); 
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/*  initialize  image  buffer  */ 

BufferlnitCkemel.indez.  keznel.aize,  image.col,  imagebuffer,  in); 

/*  display  data  in  image  buffer  before  shifting  */ 
for(i  =  0;  i  <  kemel.aize;  i++  ) 
for  Cj  =0;  j  <  image.col;  j  ++  ) 

printf  ("Image  Data  [Xd]  [Xd]  =  Xd\n",  i,  j,  imagebuffer [i]  [j]) ; 

/*  shift  1  row  of  data  in  image  buffer  */ 

BufferShlfter(kernel_lndex,  kernel.size.  image.rov,  image.col.  10,  imagebuffer,  in); 

/*  display  data  in  image  buffer  after  shifting  •/ 
ford  =  0;  i  <  kemel.size;  i++  ) 
for  (j  =0;  j  <  image.col;  j  ++  ) 

printf  ("Image  Data  CXd]  CXd]  =  Xd\n",  i,  j,  imagebuffer  Ci]  [j]) ; 
free(imagebuffer) ; 
return  0; 

} 


2.3  Image  Convolution 

DECRIPTIONs  perform  2D  convolution. 


USAGE? 

float  ConvolutionUnit(kcrnJn<lcx,  coUmage,  colJoop,  imagcJ)uf,  kcrn.buf) 


int  kem  Jndex:  index  of  Gabor  type  kernel, 
int  coIJmagc:  column  of  image. 

int  colJoop;  actual  column  number  in  convolution  double  loop, 
unsigned  char**  imageJmf:  pointer  of  a  block  storing  image, 
unsigned  char**  kcrn.buf:  pointer  of  a  block  storing  kernel. 


REMARKS: 

ConvoiutionUnit()  performs  2D  commlution  with  image  data  from  im¬ 
age  buffer  and  kernel  <lata  from  kernel  buffer.  The  column  boundary 
in  image  plane  has  been  taken  care  of  in  this  mo<lulc,  while  the  row 
boun<lary  in  image  plane  has  been  taken  care  of  by  DufferShifter(). 
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The  coUoop  U  the  actual  column  number  in  convolution  <louhle  loop. 
After  convolution  operation,  the  image  buffer  is  up<late<l  by  Buffer- 
Shifter(]  so  that  next  convolution  operation  will  be  performed  accord¬ 
ingly  with  the  input  of  right  set  of  image  data.  The  image  format  is 
8-bit  per  pixel  in  binary,  an<i  the  image  buffer  is  static  un»gncd  char 
<loublc  pointer  type. 


EXAMPLE: 


#lnclude<8tdlib .h> 

«lnclude<stdio  .h> 
tine  luda  <matli .  li> 

*lnclud«<malloc .h> 

int  mainO 

{ 

static  unsignad  char  *aimagabuffar; 

static  float  aakarii.t,  **kani_x,  aakam.y; 

float  spatial.fraq  =  .1,  oriantation  =  0; 

int  filtax.typa  =  0; 

int  row,  col,  kamsiza,  index ; 

int  inaga_rav=50,  imaga.colsSO; 

float  imaga; 

FILE  ain,  aout; 

/a  calculata  kamal  siza  and  indax  a/ 
karnsiza  =  (int)cail(1.0/(8patial_fraq)); 

if  (f  abs  (  (float )  karnsiza/2 . 0-cail  (  (float )  k9m8izo/2 . 0}  )  <=0 , 01 ) 
kamsizaaa ; 

indaxs(int)  (((float  )kamsiza-l  .0)/2.0) ; 


/a  allocate  memory  for  kernels  a/ 
karn_t  =(float  **)  mallocCkemsizeazizeof (floata)) ; 
karn.x  =(float  **)  mallocfkernsizaasizeof (floata)) ; 
karn_y  =(float  **)  malloc(kernslzaasizaof (floata)) ; 
for  (i=0:  i<  kamsiza;  i++  )  { 

kern_t[i]  =  (float  a)nkalloc(kemsizea8isaof (float}+l) ; 
karn_x[i}  =  (float  a)aalloc(kemaiseasiseof  (float)al) ; 
kem.yCi]  =  (float  a)malloc(kemBizoa8izaof  (float)-*-!) ; 
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} 

/*  allocate  namory  for  Image  buffer  */ 

imagebuffer  ^(uiialgned  char  *«)  mallocCkemalzeeslzeof (unalgned  char*)); 
for  (i=0;  i<  image.col;  1++  )  { 

Imagebuffer [i]  (unsigned  char  *)malloc(imaga_col*slzeof (unsigned  char)); 

} 

/*  generate  Gabor  type  kernels  */ 

CreateGabor(8patial_freq.  orientation,  filter.type.  index,  kem.t,  kem.x,  kern_y>; 


/*  open  image  file  in  REiD.QNLY  and  BINARY  mode  */ 
if  ((in  =  fopenC'image.dat".  ”rb"))  ==  NULL)  exit(l); 
if((out  =  fopen("conv.dat'',  "»“))  »  NULL)  exit(l); 

/*  initialize  image  buffer  */ 

Bufferlnit(kernel_index,  kemel.size,  Image.col,  imagebuffer,  in); 

/*  do  convolution  •/ 

for(  row  =  0;  row  <  image.rov;  row  ++  )  { 
for  (  col  ■  0;  col  <  image.col;  col  ++  )  { 
image  =  ConvolutlonUnlt (index,  image.col,  col,  imagebuffer,  kern.t); 
fprlntf(out,  "Xf\t",  image): 

/*  shift  data  in  image  buffer  */ 

BufferShifterdndex,  kemslze,  Image.rov,  image.col,  row,  imagebuffer,  in); 

) 

f prlntf  (out ,  ''\n") ; 

> 


free  (kem.t): 
free(kem_x) ; 
free(kem.y); 
free (imagebuffer) ; 
return  0; 

> 
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3  MOTION  DETECTION  ON  2D  TEST 
PATTERN  IMAGES 

This  section  <lcsc.ril>cs  &  motion  <ictcction  with  ono  orientation  and  spatial  fre¬ 
quency  tunc<l  Gahor  f><nction.  Unlike  the  algorithm  to  he  (Ic8crihe<l  in  the  next 
section,  this  motion  detection  is  majorly  usc<l  as  testing  purpose  for  the  know 
sine  or  cosine  wave  image  patterns,  whicli  can  he  tuned  to  only  one  orientation 
and  spatial  frequency. 

DESCRIPTION:  motion  detection  on  a  cosin  wave  moving  images. 
USAGE:  GaliorCos 


input: 

(1)  Modulated  Spatial  frequency  (wo)  of  the  Kemd  (Cycles  per  pixel) 

(2)  Length  of  the  kernel  in  sigma  value  ( length  =  l/{8igma*omcga)) 

(3)  Kernel  Orientation 

(4)  ThresoUl  for  Motion  Detection  (to  eliminate  extreme  high  value 
of  motion  (u,v)  from  the  homogeneous  image  region  where  derivative 
features  are  poorly  define<l 

(0)  Kernel  Phase  (0:  Gabor  Real  Part;  90;  Gahor  Imaginary  Port) 

(6)  Derivative  of  Image  File:  file  name  of  the  derivative  (finite  differ¬ 
ence)  of  the  image  frame  3  minus  frame  I 

(7)  Image  File:  file  name  of  the  2nd  image  frame 

output: 

(1)  u.ilat  motion  8pcc<i  along  X  axis 

(2)  v.<lat  motion  spee<l  along  V  axis 

(3)  tl.<lat;  gahor  response  to  <lerivative  of  moving  image  (Qt) 

(4)  xl.dat:  gradient  gahor  response  to  moving  image  (Dt) 

(5)  yl.dat:  gradient  gahor  response  to  moving  imiigo  (Qy) 

The  kernel  size  an<l  image  size  arc  definc<l  in  #dcfine  macros. 
REMARKS: 


This  program  extracts  an  motion  vector  hy  Grwlient  Gal)or  Model.  A 
set  of  three  consecutive  image,  frames  is  nccdnl  for  extracting  motion 
information.  The  first  fname  nn<l  the  third  frame  arc  usc<l  to  form 
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imago  derivative  l)y  finite  <liffcrcncc  tediniquc,  then  the  (liffcrcncc  is 
convolvc<l  with  the  Gal>or  kernci  to  produce  Qf  The  second  frame 
is  convolved  with  derivatives  of  Gai>or  (in  x  nn<l  in  y  direction)  to 
produce  Q,  and  Hy.  The  user  lias  to  first  define  the  orientation  and 
spatial  frcquoncy  of  the  cosine  test  pattern.  The  selection  of  the  length 
of  the  Gahor  kernel  has  to  he  made  with  the  consiilcration  to  the 
orientation  and  spatial  frequency  selectivity.  As  a  general  guideline, 
the  user  is  a<lvisc<l  to  use  the  following  formula  to  estimate  tlie  Icngt.h 
value; 

Length  -  ^ 


where  I!  is  the  sigma  of  the  test  pattern. 

Another  point  for  the  user  to  realize  is  that  tlircsholding  is  performed 
on  the  processe<l  image  ilata,  Gc,  Gy,  and  Qt,  to  deal  with  very  small 
value  at  certain  pixel  location.  With  this  thrcshohling  the  calculatcfl 
optical  flow  vectors  at  these  locations  vrill  he  free  of  erroneous  results 
due  to  the  nature  of  the  early  virion  computing  (  most  of  them  arc 
mathematically  characterizeil  as  ULposeil  problems). 

MATHMATICAL  FORMULATION: 

n*  =  n,u + Oyw 

where  is  the  convolution  of  Gahor  kernel  with  derivative  of  image,  R.  is  the 
convolution  of  x  gradient  Gahor  kernel  with  moving  image,  Ry  is  the  convolu¬ 
tion  of  y  gradient  Gahor  kernel  with  moving  image. 


EMAMPLE:  CosConv 

Modulate<l  Spatial  frequency  In  KcmcI(Cyclcs/pu:cl):  0.1 
Length;  4.0 

Kernel  Orientation(<lcgrce):  30.0 
Thrcsohl  of  Motion  Detection:  100.0 
Kernel  Phasc(dcgrec);  90.0 
Derivative  of  Image  File:  <lcvi.img 
Moving  Image  File:  movc.img 
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4  MOTION  DETECTION  ON  2D  IMAGES 
WITH  MULTIPLE  KERNELS 

Far  any  realistic  motion  detection,  different  spatial  frequency  and  orientation 
features  shoul<l  all  be  considered.  This  requires  the  use  of  more  than  one  spatial 
frequency  an<l  orientation  tuncel  Gahor  kernels.  This  section  describes  how  this 
multiple-kernel  process  can  be  accomplished. 

DECRIPTION:  Motion  Detection  on  Moving  Images. 


iincluds  "Isaatst.e" 


USAGE:  GaborMotion 


input: 

(1)  Morlulaterl  spatial  frequency  (wo)  of  the  kernel  (cycles  per  pbtel) 

(2)  Starting  orientation  of  the  komel  (degree) 

(3)  Rotating  angles  each  time  (degree) 

(4)  Number  of  rotations,  or  the  number  of  different  oriented  Gabor 
kernels 

(5)  Thresohi  of  motion  detection  (a  rccommcndated  value  is  greater 
than  00) 

(6)  Kernel  phase  (0:  Gal>or  R.cal  Part;  90:  Gnl>or  Imaginary  Part) 

(7)  Derivative  of  Image  File:  file  name  of  derivative  (finite  difference) 
of  the  average  of  image  frame  3  minus  frame  1  (8)  Moving  Image  File: 
file  name  of  image  frame  2 

output: 

(1)  u.dat  motion  speed  along  X  axis 

(2)  v.dat  motion  spee<l  along  Y  axis 

The  kernel  size  an<l  image  size  arc  defined  in  #dciine  macros. 
REMARKS: 


This  program  cxt.raet.s  an  optical  flow  from  3  consecutive  image  frames 
by  Gradient  Gabor  Model.  Wlien  selecting  iliffcrcnt  length  of  kernel. 
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the  con6i<icrat.ion  has  to  he  given  in  the  connection  to  ttic  preprocess¬ 
ing  of  l)an(lpnss  filtering  or  lowpass  filtering  operation.  We  suggest 
the  use  of  Gaussian  kernel  for  lowpass  filtering  or  LoG  (Laplacian  of 
Gaussian)  kernel  for  l)an<lpass  filtering  before  the  use  of  Gabor  kernels 
for  motion  <lctec.tion.  As  one  general  guiilclinc,  the  same  sigma  value 
of  the  Gabor  kernel  can  be  used  for  the  Gaussian  kernel,  and  the  sa^nc 
sigma  of  the  Gabor  kernel  can  be  used  for  the  LoG  kernel. 

MATHEMARTCAL  FORMULATION: 

+  n*1i  -I-  QyV  =  0 

where  Qt  =  . ft,  =  [fti,ft2,...,ft‘f ,  ft^  =  (ftj,ftj . ftjf , 

n{=  is  the  convolution  of  Gabor  kernel  with  <lcrivative  of  image, 

ftj  is  the  convolution  of  x  gra<lient  Gabor  kernel  with  moving  image,  and 

fty  is  the  convolution  of  y  gradient  Gabor  kernel  with  moving  image. 

EXAMPLE:  GaborMotion 


Modulated  Spatial  frequency  In  Kemcl(Cycles/pixel):  0.1 

Stating  Orientation  of  Kerncl(<ltgrcc):  0.0 

Rotating  Angies  Eacli  Time(<legrcc):  30.0 

Number  of  R.otations:  4 

Thrcsold  of  Motion  Detection:  100.0 

Kernel  Phasc(dcgrec):  0.0 

Derivative  of  Image  File:  dcvi.img 

Moving  Image  File:  movc.img 


5  LEAST  SQUARE  ESTIMATION 


DES  CRIPTION:  Least-Square-Estimation 

USAGE:  solve2(int  row,  int  column,  float  thresohl) 
or  solve4(int  row,  int  column,  float  thrcsold) 


int  row:  the  row  of  image; 

int  column:  the  column  of  the  image; 

float  thrcsold:  the  thrcsold  of  motion  detection. 
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REMARKS: 


This  program  performs  a  toasf.-squarc-csMmation  to  minimize  the  q- 
ror  of  optical  flow  computation.  The  tediniquc  will  allow  the  better 
han<lling  of  the  ili-po8e<l  problem,  optical  flow  computation.  As  it  is 
well  known,  many  early  vision  problems  arc  ill-pose<l  in  nature.  To  be 
able  to  (leal  v  th  this,  a  least  square  estimation  has  been  developed  for 
optical  flow  computation  with  a  set  of  2  elifferent  orientation  Gabor 
kernels  or  a  set  of  4  different  orientation  Gabor  kernels  respectively. 


MATHMETICAL  FORMULATION; 

E(u, «)  =  +  fi"»u  +  O’yuj 

EXAMPLE;  solve2{60,  60,  200.0) 


SD  Motion  Detection 


25 


The  cosine  pattern  was  moved  .2  pixels  at  30.5  degree  orientation 
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0  5  10 

X  Axis  VI 

Average  value  of  the  computed  move:  1.97;  Imaginary  Gabor  kernel 


Figiirn  3.  Mntimi  Hfitnction  rmiilt  for  roniiui  wovo  iinngon,  thrm 
fromiw  went  geiioroted  with  nrionUtion  30.5  dogrm  mid  npniiiil  froqiii-nry  0.025. 
The  motion  detor.tion  wn*  performed  by  Gebor  kernel  with  orientetion  luid 
epatiel  fre<|ui«r.y  equal  to  30.5  degree  and  0.025  ranpeetively. 
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THREE-DIMENSIONAL  MOTION 

DETECTION 


ABSTRACT 

1  VIRTUAL  REALITY  TECHNIQUE  FOR 
CREATING  3D  TEST  PATTERNS 

To  he  able  to  test  the  motion  detection  algorittim  icalu^tically  with  better  con¬ 
trol  of  the  testing  comlition,  we  have  dcvdopeil  the  tccliniquc  based  on  the 
concept  of  virtual  reality  to  creat  virtual  camera  model.  We  have  been  using 
the  tool  developed  by  Watkins  A  C-type  language  was  used  to  creat  the 
user  desirc<l  testing  conditions  wliicli  include  the  user  selection  of  the  following 
features; 


■  A  world-coonlinate  system  where  every  objects  arc  defined  and  arc  refer¬ 
enced  to,  an«l  a  viewcr-coor<linatc  system  v-liidi  is  observer-centered  coor¬ 
dinate  system,  as  well  as  the  relative  position  between  the  two. 

■  A  camera  mo<lel  with  the  choice  of  orientation  (pointing  direction)  of  the 
camera,  the  optical  characteristics  of  the  lens. 

■  Lighting  con<lition  with  user  control  of  the  placement  and  the  number  of 
point  light  sources. 

■  3D  object  construction  using  primitive  building  blocks,  and  the  clioice  of 
<lifferent  reflection  models,  such  as  <lifFuscd  reflection,  specular  reflection, 
and  ambient  light, 

C.D.  Watkins,  S.D.  Coy,  aiiH  M.  Finlay.  PhotorrAli/m  onrf  liny  Tracing  in  C,  M&T 

Books,  Now  York,  NY  10011,  1992. 
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Given  in  tlic  floppy  disk  numi)cr  2,  in  directory  MOVIE,  one  of  tiic  image 
sequence  gcnerate<i  by  this  tedmique,  ONR.1.FLI,  can  be  displayeil  as  a  movie 
on  486/50MHz  machine  witii  SuperVGA  monitor.  To  <lisplay  this  sequence, 
first  load  t)ic  floppy  disk  into  a  floppy  disk  <iriver,  tlien  go  to  this  drive  an<i 
go  to  directory  MOVIE,  tiien  type  PLAY  ONIl,l.FLI  to  execute  tiie  display  of 
the  movie.  Please  be  aware  of  the  software  copyriglit  issue,  the  program  in  this 
floppy  disk  shouhl  not  be  copyed,  or  used  for  any  other  purpose.  To  obtain 
a  copy  of  this  utility  program,  please  contact  the  publisher  indicated  by  the 
program. 

With  this  virtual  reality  technique,  the  user  can  creat  well-controlled  motion 
images  with  the  exact  known  camera  setup  and  the  relative  position  between 
the  observer  and  the  moving  ob  jects.  This  is  a  very  powerful  way  of  testing 
and  verifying  the  algorithm  <levcloped  for  motion  computation  and  motion 
compensation. 


2  CREATING  IMAGE  SEQUENCE  WITH 
BOTH  OBSERVER  AND  OBJECT 
MOTION 

Using  the  tccliniquc  <lcscribcd  in  the  previous  section,  the  user  can  creat  image 
sequence  with  both  ohserver  and  object  motion.  This  section  describes  one 
example. 

DESCRIPTION:  generate  squence  of  moving  images  including 


(1)  sphere  moves  along  Z  asix  with  step  2  units; 

(2)  observer  moves  along  Y  axis  with  step  2  units; 

(3)  combination  of  (1)  and  (2). 


The  program,  BMll.B,  for  creating  3D  virtual  environment  anil  the  image 
sequence  is  given  the  last  chapter  of  this  manual. 
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3  3D  MOTION  DETECTION 

Using  tlic  program  GAR.DOMOTION.  tlirco  dimensional  nioi.ion  can  be  com¬ 
puted.  Tile  prngi-am  is  i.lic  int.egrat.ion  of  several  program  modules  dcscrilml  in 
Chapter  3.  A  source  code  of  the  program  is  given  in  Clia|ii.ei  5.  The  program 
wat.  documented.  With  the  understanding  of  tlie  programs  in  Chapter  3,  the 
use  of  this  piogi  am  shall  not  pose  new  challenge. 


Enlarged  portion  without  showing  the  x-y-z 
coordinate  (movie  file  is  included  in  the 
floppy  disk  and  can  be  played  on  486/50Mhz 
machine  or  above  with  SVGA) 


4.  .A  fuitiir  t*f  ;i  *•  of  iinaj'o.'*  wii-li  i'l»!«'*rver  Jind 

nphor^*  ;uul  w>  i<'  iti'itcribAd 

in  t-ln*  .i'  dirr..  I  <  III. •  M  1  1 . D 


Z  Direction 
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Y  Direction  ( Observer  Moves  5  Units  Along  Y  Axis) 
30|  '  '  "I  I  I  I 


20 


10 


S 


0 


-10 


-20 


_ \ - 1 - 1 _ I _ u. 

-20  -10  0  10  20 
Kernel  31x31 , 4  Orientations,  35  Degrees  Apart.  Phase  0 


Figure  S.  Motion  detection  reniilt  for  a  »phere  image  with  4  kernels 
45  degree  apart. 
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ABSTRACT 

This  cliaptec  gives  tltc  source  co<le  listing  of  all  the  programs  ilescriiieil  in  this  manual. 
A  MS-DOS  floopy  disk  is  also  provided. 


1  BARIMAGE.C 

/♦ 

a  Copyright  (c)  Dr.  Hua  Li’s  Lab  1994.  All  Rifats  reaorved. 
a  Purpose;  To  generate  an  oriented  bar  stripe  laage  pattern  for  verlflng 
a  the  concept  of  spatial  frequency  and  orientation  bandvldth 

a  selectlvitias. 

a  Image  Format:  binary  image  from  0  to  3S5.  The  sine  is  ROU  z  COLUMN; 
a  Name:  Bar Image. c. 

a  Compilation:  (g}cc  -o  outputfllenaae  Barlmage.c  -Im 
a  Side  Effect;  No. 

a  Input;  Spatial  frequency.  Orientation,  and  slsa  of  image 
a  Output;  bar  image  file, 
a  Coded  by:  Xiaohul  Heng, 
a  Executed  Machine:  Sun  Workstation 
a  Last  Change:  Feb.  24,  1994 
a/ 


#lnclude<strlng . h> 
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#liicludo<8tdio  .h> 

#iiiclud«<nath  .li> 

Barlmgdnt  row,  Int  colium,  float  omega,  float  orient,  char*  filename) 

•C 

FILE  efile; 
int  i,  j,  image; 

float  pi=3.1416,  Spatial_Freq_X,  Spatial_Freq_Y,  Orientation,  value; 

Spatial_Fxeq_X  =  omega  *  coaCorient  *  pi  /  180.0); 

Spatial.Freq.Y  &  omega  *  sinCorient  e  pi  /  180.0); 

if((file=fopen(filename,"wb"))==NULL  )  { 
printfC"  Can’t  open  output  file  for  writing"); 
exit(-l) ; 

} 


/* 

*  The  bar  wave  image  formulation: 

*  value  =*  255  if  cos(2*pi*(x*wx  +  y*wy))  >0.0 

e  0  otherwise 

*  where  wx  is  the  spatial  frequency  component  in  u 

«  wy  is  the  spatial  frequency  component  in  v 

*  then  the  value  of  image  is  mapped  to  0-265 

•/ 

for  (  i  =  0;  i  <  row;  1++)  < 
for  (  j  =  0;  j  <  column;  j++)  { 

if{  cos(2.0*pie(i*Spatial_Freq_X  +  j*Spatial_Freq_Y))  >  0.0  ) 
image  s  265; 

else 

image  =  0; 

fprlntf  (file ,  ”Xd\t''  .image) ; 

> 


fprintf(file,"\n"); 

> 

fclose(file) ; 

} 
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2  COSIMAGE.C 


Copyright  (e)  Dr.  Hua  Li's  Lab  1994.  ill  Rights  Rassrvad. 

Purpose  To  gsnarata  a  orientod  cosia  moving  image  and  image  derivative. 
This  pattern  is  used  to  test  the  coordinate  system  rotation 
and  image  flov. 

Name  :  Coslmage .  c . 

Compilation:  cc  -o  outfilename  Coslmage. c  -Im 
Output  files:  (1)  orig.img  —  original  image. 

(2)  move.lmg  —  moved  image. 

(3)  devi.img  —  derivative  of  image. 

Coded  by:  Xiaohui  Meng. 

Last  Change:  Peb.  23  1994 


tinclude  <stdio.h> 

*include  <stdlib.h> 
iinclude  <math.h> 

CosImgCint  row,  Int  column,  float  omega,  float  orient,  float  speed,  float  time) 

{ 

FILE  *infilel,elnfile2,elnfile3: 
float  value 1, value 2, values; 
float  vz,vy,pl=3.1416; 
int  i,j; 


/•  Image  Pattern  Generator 

*  The  formula  of  cosin  image  la: 

♦  I(x,y)*  co8(2epi*omegae(x*cos(thata)'*'yesin(theta)-8peed*tlme)) . 

•  where  omega  is  the  spatial  frequency 

*  theta  is  the  orientation 

e  the  product  of  time  and  speed  defines  the  phase  offset  along  theta 

•  direction. 

*/ 

vx=cos(theta*pi/ 180.0) ; 
vy=sln(thetaepi/180.0) ; 


if (Cinfilel=fopen("orig.img","wb"))==NULL)  exit(-l) ; 
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if  ( (inf  il«2=f  opeii("d«vi .  iag" ,  ''wb'OJsamJLL)  txit  (-1) ; 
if((infil«3=fop*n("m6v*.iag“,»wb"))==NULL>  •xit(-l>: 

for  (  i  =-(row-l)/2;  i<(ro»-l>/2+l;  i++) 

{ 

for  (ja-(colUBm-l)/2; j<(eoliuBn-l)/2+l; j++) 

{ 

valuol  =(int)2SS.0*(eo8(2.0*pi*om*gn*(vx*(float)i+vy*(flont)j-8p*«d*tiiM))+l)/2.0; 
valua2  =(int)26S .0*(co8(2. O*pi*oaega*(vx* (float) i+vy*(float)j-apaad*(tima+l)))+l)/2.0 
valuaS  =(int)265.0*(co8(2.0*pi*omaga*(vx*(float)i+vy*(float)j-8paad*(tina+2)))+l)/2.0 
fprintf (infilal,“Xd\t".valual) ; 
f pr intf (inf ila2 , "Xd\t" , (valua3-valua 1 > /2 . 0) ; 
fprintf  (infil«3.'*Xd\t'*,valua2} ; 

} 

fprintf (inlilai ."Vn") ; 
fprintf(infila2,"\n"): 
fprintf  (inf  ile3,"\n’') ; 

} 

fclosa(infilal) ; 
fclo8a(infila2) ; 
fclosa (infilaS) ; 

} 


3  COSCONV.C 


/• 

*  Copyright  4  Dr.  Hua  Li's  Lab  1994.  All  Rights  rassrvad. 

*  Purposa:  To  calculata  tha  imaga  flow  from  moving  oriantad  cosin  imagas. 

a  Kamal  :  Gabor  raal  part  karaal  and  inare  sisa  ROWSIZE  x  COLUMNSIZE. 

*  Imaga  Format:  Binary  Fila 

a  Nama  :  CosConv.c. 

a  Compilation;  cc  -o  outputfilanama  CosConv.c  -Im 

a  Slda  Effact;  Llmltad  by  maximum  array  sisa 
a  Input:  (1)  Xarnal  Modulation  Spatial  frocpiancy 
a  (2)  langth  (  langth  s  l/(slgaaaomaga)) 

a  (3)  Karnal  Oriantatlon 

a  (4)  Thrasold  of  Motion  Dataction 
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(5)  K<ro«l  Phas*  (0:  Gabor  Real  Part:  90:  Gabor  Imaginary  Part) 

(6)  davi.img:  fila  nano  of  dorlvatlva  of  moving  Image 

(7)  mova.img:  file  name  of  moving  Imago 
Output: (1)  u.dat  motion  apood  along  X  axla 

(2)  v.dat  motion  apood  along  Y  axla 

(3)  tl.dat:  gabor  roaponao  to  dorlvatlvo  of  moving  Imago 

(4)  xl.dat;  gradlont  gabor  roaponao  to  moving  Imago 

(5)  yl.dat:  gradient  gabor  roaponao  to  moving  imago 
Exocutod  Machine:  Sun  Sparc  Station 

Coded  by:  Xiaohul  Kong. 

Laat  Change:  Kay.  25  1994 


•Include  <8tdio.h> 
•include  <atring.h> 
•include  <8tdllb.h> 
•include  <math.h> 


•define 

KERNELSIZE 

El 

/*  reaerved  kernel  aixe  */ 

•define 

ROWSIZE 

81 

/* 

roH  aixe  of  image  */ 

•define 

COLUKNSIZE 

81 

/* 

column  alze  of  image  «/ 

void  mainO 

{ 


float  rparttCKERNELSIZE]  [KERNBLSIZE]  ,rpartx[KERNELSIZE]  [KERNELSIZE] ; 
float  rparty [KERNELSIZE]  [KERNELSIZE], 

bufforl [KERNELSIZE] [COLUHNSIZE] ,buffer2 [KERNELSIZE] [COLUKNSIZE] ; 
FILE  *lnfllol,oinfile2; 

FILE  ooutl ,*out2,eout3,«out4,*out5: 

Int  row, col .rovl, coll ,rov2,col2,l,j , index, nufflbor,komslza; 
f 1 oat  vx , vy , omega , theta, pl=3. 1416, algma , orient , length , thr eaold ; 
float  image , image 1 , lmage2 , images .phaae ,u , v ,urof , vref ; 
char  dovi.name  [20] ,  move.name  [20]  ; 

printf("\n  Input  Spatial  FrequencyCCyclaa/plxel) :"} ; 

acanf ("Xf ",Romega) ; 

printf("\n  Input  length 

acanf  ("Xf  "  ,4Elength) ; 

printf("\n  Input  Orientation 

acanf ("Xf" ,ktheta) ; 
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prlntf("\n  Input  thrtnold 
Bcanf  ("Xf  "  .Ictliraaold) ; 

pxlntf("\n  Input  Phns*  alilftlng(d«gr*«B)  ; 
scan!  ("X<"  .Jtphas*) ; 

prlntf("\n  Input  Darlvntiva  Imng*  Flla:"); 
scanf  C'Xs" .  dflvl.nama) ; 
ptintf("\n  Input  Moving  Imag*  File;"); 
ncanfC'Xs",  aove_na0e); 

sigma  =  1.0/(?.ength*omega) ; 
kemslze  =  (int)ceil(l  .0/(omaga)> ; 

if (f abs ( (float }kernsize/2 . 0-ceil ( (f loatlkazusiza/S . 0) ) <-0 .01) 
kentsize+'t'  • 

pzlntf ("KernSlzesXd\n" .kernsize) ; 
indez«i<int)C((float)kernslza-1.0)/2,0) ; 


if ((infilel®fopea(devi_aame,"r"))**NULL)  exlt(-l) ; 
if  ((infile2=fopen(move_naiBe,"r"))«NULL)  exit(-l) ; 


/* -  Gabor  type  kernel  generator  - - - */ 

vxsonegaecos(thetaepl/180) ; 

vysomega^ainCthetaapi/lSO) ; 
for  <i=-index;i<index+l;i++)  { 
for  (js-ittdex;j<lndex+l; j++)  { 

imagesexp(-(float)(i«l-*'j*j)/(sigma*8igBa)) ; 

imagelscoa(2.0*pie(vxe(float)i-^vye(float)J)-«'phase*pl/lB0.0)*image; 
iinage2ssin(2 .  OapiaCvz^Cf  loatli-t-vyeCf  loat)  J  )'4’pba8e*pl/160 ,0)*image ; 
rpartt  [i-t-index]  [j-^indezl^iBagel; 

rpartzCl-t-indez]  [j'»-index]=-iBage2*Txv2.0«pl-iBagala(float)l*2 .0/(aigmaeaigna) ; 
rparty  [i-t-indez]  [J  -t-index]  =-image2*vy«2 .  Oepi-iaaga  1*  (f  lost )  j  *2 . 0/  (aigmaeaipaa) ; 

} 

> 

/♦ - convolution - */ 

if ((outl=fopen("tl .dat",*V))==NULL)  exit(O) ; 
if ((ottt2=fopen("xl.dat","w"))==NULL)  exit(O) ; 
if ((out3=fopen("yl.dat","w"))=sNULL)  exit(O) ; 
if ((out4=fopen("u.dat","w"))==HULL)  oxit(O) ; 
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if ((outB»fop«n("v.daf.«w«))»«MOLL)  •xit(O) ; 

fox  (rovBiad«x;rov<k«ra«ls«;Tov^)  /*  Inltlxlixa  tho  buff at  */ 

{ 

for  (co1sO:co1<CQLUKNSIZE:co1+4-) 

f scar  f  (inf ila  1 ,  "Xf '*  ,kbuf f  arl  [row]  [col] ) ; 
facanf (inflla2, "Xf " ,kbuff ar2[rott] [col] ) ; 
bttf  f  arl  [rov-indax]  [collcbuffarl  [indax]  [col] ; 
buf f ar2 [row-indax] [col] =buff 0x2 [indax] [col] ; 

> 

} 

for (i=0;i<R0USIZE;  /*  do  tba  convolution  */ 

i 

for  (  j=«0;j<C0HmMSIZE;j++) 

< 

imago laO.O; 
lmaga2=0.0; 
imaga3«0 . 0 ; 


for  (rowl*-indax;rowl<indax+l;rowl++)  /•  convolution  unit  ♦/ 

■i 

for  (colls-indax;coli<indox4-l;coll-*"*-) 

{ 

col2=j+coll; 
if(col2<0)  col2=0; 

if  (col2>C0LUMNSIZE'l)  col2=C0LimNSIZE-i ; 

imagal-)-=buffarl[ro«l>indox]  [col2] arpartt  [indax-t-rovl]  [indax^coll]  ; 
imaga2-)'=buf f 0x2 [rowl-»^indax]  [col2]«rpartx[lndox^xowl]  [indox+coll]  ; 
imagaS-t-sbuf f  0x2  [rovl-'-indax]  [col2]arpBrty[indox-*-xovl]  [indox'^coll]  ; 

> 

} 


/• 

a  do  motion  spaad  calculation 

•/ 


if(  faba(imaga2)  >  tbxoaold  ) 

u  ~  cos(pi*tbata/180.0)aco8(piatbata/180 .0)'*l]Bagal/lmaga2; 
<:l«a  u  =  0.0; 
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lf(  fabs (images)  >  thresold  ) 

V  =  sin(pi*tbeta/180.0)*5ln(pi*theta/180.0)*linagel/image3; 
else  V  =  0.0; 

/♦  if(uref  =0.0)  u  =  0.0; 

else  u  =  uref*(vref*vref/(uref*uref+vref*vref )) ; 
ifCvref  =  0.0  )  v  =  0.0; 

else  V  =  vref*(uref*uref/(uref*uref+vref*vref)) ;  */ 

fprintf (outl .image 1) ; 

fprintf (out2 ,"Xf \t” ,image2) ; 

fprintf (outs, images) ; 

fprintf (6ut4 ,"Xf \t" ,u) ; 

fprintf  (outs  ,"j4f  \t" ,  v)  ; 

}  /♦  for  loop  for  j  */ 

fprintf  (out  1 , \n" ) ; 
fprintf (out  2 , " \n " ) ; 
fprintf (outs," \n") ; 

fprintf  (out4,"\n*‘) ; 
fprintf  (outs, ''\n") ; 

for  (col=0;col<C0L0MNSIZE;col++)  /*  shifter  •/ 
for  (rou=0;row<kemsize-l;row++) 

•C 

bufferl  [row]  [col]=bufferl[row+l]  [col]  ; 
buffer2[rov]  [col]=buffer2[row+i]  [col] ; 

} 

if  (i<R0WSIZE-index-l) 

< 

fscanf  (infilel,"y,f",*bufferl[kem8ize-l]  [col]) ; 
f8canf(infile2,"Xf",kbuffer2[karn8lze-l]  [col]) ; 

> 

> 

} 

f close (outl) ; 
f close (out2) ; 
f close (outs) ; 

fclose(out4) ; 
fclo8e(out5) ; 
fclose(infllel) ; 
i(-lose(infile2)  ; 
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> 


4  CREATGABOR.C 

/* 

*  Copy  right  Dr.  Li’s  Laab  1994.  All  Rights  rassrved. 

*  Purposs:  To  gensrata  Gabor  kemal. 

*  Format :  float  point 

*  Name:  CraatsGabor .c 

*/ 

void  CreatGabor(omaga,  thata,  typa.size,  tampt,  tampx,  tempy) 
float  omaga,  theta; 

Int  typa,  slza; 

float  **teiipt,  •♦tampx,  **tampy: 

{ 

float  omaga.x,  omaga.y,  sigma; 

float  phasa  =0.0;  /*  Gabor  raal  part  as  default  */ 

float  image,  imagal,  image2; 
float  pi  =  3.141593; 
int  i,  j; 

omaga_x  =  omega*cOB(thataapi/180) ; 
omaga.y  =  omaga*8lu(thata*pl/180) ; 
if (typa  ==  1  )  phasa  =  90.0; 
else  phasa  =  0.0; 

8lgma=  1 . 0/ (4 . 0*omega) ; 

for  (i=8ize ;i>=-8iza ;1 — ) 

for  (j=8iza; j>=-size; j — ) 

■C 

imaga=exp(-(float) (i*i+j*j)/(0igma*8igma)) ; 

imagal=cos  (2. 0*pi-<' (omaga.x*  (float>i+omega..ya  (float  )j)+pha8e*pi/180.0)*imaga; 
imaga2=sln (2. 0*pl* (omaga.x* (float)i+omegB_y*(float}J)*pha8e*pi/180.0)*lmaga; 
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t«Bq>t[8iza~l]  C8lz«-j]=lfflag«l; 

teinpx[slz«-i]  [slza-j]  =-lBaga2*ODega_x4'2.0«pi-lnagel* (float )i*2.0/(slgiiia4'aigma) ; 
tainpy[size-l]  Cslze-j]=-iDaga24'Oia«ga_y*2.0*pi-lnag«lo(float)J*2.0/(8igma*sigBa) ; 
> 

} 

} 


5  BUFFERINIT.C 


/* 

* 

* 


* 

*/ 


Copy  Right  A  Dr.  Li '8  Lib.  All  Rights  R«Btrv«d 
Purpo8«:  initiallz«  th«  iaage  buffer  to  take  care  of 

image  boundary  condition  for  convolution. 
Module  Name:  Bufferlnlt.c 

Required;  file  vae  opened  in  "RB"  mode. 

Image  Format:  8-blt/plxel  binary  format. 

Coded  by;  Xlaohul  Meng 

Last  Change:  Oct.  1994. 


the 


void  Bufferinit (index,  kemslze,  column,  buffer,  fp) 

Int  index,  kemaize,  column; 
unsigned  char  «*buffar; 

FILE  efp; 

{ 

int  rou ; 

for(  row  =  0;  row  <  kemsize;  row++)  { 
if(  row  <  index  )  •( 

freed (buffer [row] ,sizeof (unsigned  char),  column,  fp); 
rewind (fp) ; 

} 

else 

fread (buffer [row] ,  sizeof (unsigned  char),  column,  fp) ; 

} 
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6  BUFFERSHIFTER.C 


h 

*  Copy  Right  4  Dr.  Li's  Lib.  All  Rights  Rsservsd 

*  Purpose:  shi^t  data  in  image  buffer  for  convolution. 

*  Module  Name:  BufferShifter .c 

*  Required;  file  was  opened  in  "RB"  mode. 

*  Image  Format:  8-bit/pirel  binary  format. 

*  Coded  by:  Xiaohui  Meng 

«  Last  Change:  Oct.  1994. 

•/ 


void  BufferShifter (kern.index,  kern.size,  row.lmg,  col.ijng,  rotf.loop,  buffer,  fp) 

int  kern.index,  kern.size; 

int  rov.img,  col.img,  row.loop; 

unsigned  char  eebuffer; 

FILE  *fp; 

{ 

int  row,  col; 

for(  col  =  0;  col  <  col.img;  col++> 

for(  row  =  0;  row  <  kem.size-i;  row++) 
buffer  [row]  [col]  =  buffer  [row+1]  [col]  ; 

if(  row.loop  <  row.img  -  kern.index  -1  ) 

fread(buffer[kern_Bizo  -1] .size of (unsigned  char),  col.img,  fp); 


> 
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7  CONVOLUTIONUNIT.C 


/* 

*  Copy  Right  t  Dr.  Li ’a  Lib.  All  Rights  Reaervod. 

*  Purposo:  convolve  image  viith  kernel. 

*  Hodule  Name:  Convolutiontlnit.c 

*  Required:  file  was  opened  in  "RB"  mode, 

e  Image  Format:  8-bit/pixel  binary  format. 

*  Coded  by:  Xiaohui  Hang 

Last  Change:  Qct.  1994. 

*/ 


float  ConvolntlonUnlt(kern_lndex,  col.imags,  col.loop,  iaiage.buf,  kem.bnf) 
int  kem.lndex,  col.lmage,  col_loop; 
unsigned  char  o^lmage.buf ,  *okem_buf: 

int  rovl,  coll,  col2: 
float  conv.datas  0.0; 

for  (rowl  ®  -kem.index;  rowl  <  kaxn.indax  +  1;  rowl++  )  •{ 
for  (coll  =  -kem_lndex:  coll  <  kem.lndex  +  1;  coll++)  ■( 
col2  =  col.loop  +  coll; 

lf(  col2  <  0  }  col2  s  0;  /e judge  boundary  in  coluian  of  image  buffer  */ 
if(  col2  >  col.image  -  1)  col2  =  col.lmage  -  1; 
conv.data  -•-=  (float) image.buf  [rovl+kem.index]  [col2]  * 
kem.buf  [kem.index+rovl]  Ckem.lndox-i’coll] ; 

> 

} 

return  conv.data; 

> 


8  GABORMOTION.C 

/* 

*  Copyright  0  Dr.  Hua  Li's  Lab  1994.  All  Rights  reserved. 
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*  Puzpoae:  To  calculat*  tha  image  flow  from  moving  images. 

*  Kernel  :  Gabor  type  kernel (real  or  imaginary  part) . 

*  Image  Format;  Binary  image  vltb  RDUSIZE  x  COLUMNSIZE. 

e  Kame  :  Gabor. c. 

*  Required:  Leastsq.c. 

*  Conciliation:  cc  -o  filename  Gabor. c  -Im 

*  Side  Effect:  Limited  by  maximum  array  size. 

*  Input:  (1)  Kernel  Modulation  Spatial  frequency 

*  (2)  Starting  orientation  of  kernel 

*  (3)  Rotating  angle  eacb  time 

*  (4)  Numbers  of  rotating 

*  (5)  Threshold 

*  (6)  Phase  in  kernel.  0:  real  part;  90:  imaginary  part. 

e  (7)  Derivatevl  Image  File:  file  name  of  derivative  of  moving  image 

*  (B)  Moving  Image  File:  file  name  of  moving  image 

*  Output :(1)  u.dat  motion  speed  along  X  axis 

*  (2)  v.dat  motion  speed  along  Y  axis 

*  (3)  tl.dat:  gabor  response  to  derivative  of  moving  image 

*  (4)  xl.dat:  gradient  gabor  response  to  moving  image 

«  (5)  yl.dat:  gradient  gabor  response  to  moving  image 

*  Executed  Machine:  Sun  Sparc  Station 

*  Coded  by;  Xiaohul  Meng. 

*  Last  Change:  Mar. 16  1994 

•/ 


iinclude  <stdlo.h> 
tinclude  <string.b> 
iinclude  <8tdlib.h> 
iinclude  <math.h> 
iinclude  "Leastsq.c" 


idefine 

KERNELSIZE 

SI 

/*  reserved  kernel  size  */ 

idefine 

ROWSIZE 

61 

/*  row  size  of  image  «/ 

idefine 

COLUMNSIZE 

61 

/e  column  size  of  image  */ 

void  mainO 


float  rparttCKERNELSIZE]  [KEllNELSIZE]  ,rpartx[KERNELSIZE]  [KERNELSIZE]  ; 
float  rpartyCKERNELSIZE]  [KERNELSIZE] , 

buffer  1  [KERNELSIZE]  [COLUMNSIZE]  , buff er2 [KERNELSIZE]  [COLUMNSIZE]  ; 
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FILE 

FILE  «outl,*out2i*out3: 

int  row, col ,rowl , coll  ,row2 ,eol2,l , j , indoi .numbor .komslzo ; 

float  vx,vy,oin«ga,theta.pl=3 . 141B926B36,Bi®na,angla_start .aagla.tuna .oriaat.nutti 

float  thro sold; 

float  imago , Imago 1 , lmago2 , imago3 .phaso ; 
char  buffor_a[lE] ; 
char  buffor_bCl&] ; 
char  buffor_c[15]  ; 
char  buffor.dClB]  ; 
char  •8tringl=".dat"; 
char  *8tring2®"t\0",* 
char  ♦8tring3="x\0"; 
char  *8triiig4="y\0": 
char  dovl_iLamo[20] ,  moTO_nafflo[20]  ; 

printf("\a  laput  Spatial  FroquoncyCCyclos/pixol) ; 

8caaf ( "Xf ” ,fcomoga) ; 

priatf("\a  Input  laitial  AnglosCdogrooa) : ") > 

8canf ("Xf " .ftanglo.start) t 

printfC'Va  Input  Tuning  Anglos (dogroos) ; 

scanf  (  "Xf  •'  ,Aanglo_tuno) ; 

printf("\n  Input  Oziontatioa  Numbors 

scanf ("Xf" ,ftorlont_Bum} ; 

prlntf<"\n  Input  Throshhold  :">; 

scanf ("Xf" lAthrosold) ; 

printf("\n  Input  Phaso  shlfting(dogroOB} ; ") ; 
scanf  ("Xf'.Aphaso); 

priatf("\n  Input  DorlvatiTo  Imago  Fllo:"); 
scanf ("Xs",  dori_aamo); 
printf("\n  Input  Mowing  Imago  Filo:"); 
scanf ("Xs",  mowo.aamo); 

numbor=0;  /*  numbor  of  oriontations  */ 

sigma-  1.0/(8.0*omoga):  /*  sot  longth  ^  4.0  */ 

koxnsizo  =  (lnt)coil(1.0/(omoga}) ; 

if (f abs ( (float )kom8izo/2 . 0-coil ( (float )kotn8izo/2 . 0) ) <=0 . 01) 
komsizo-t’-t-; 

printf  ("KoxnSlzosXd\n",komsizo) ; 
indox= (int ) ( ( (f loat)koxnsizo- 1 .0} /2 . 0) ; 
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if ((infile 1=1 op«n(d«vl_n«m*,"r">)==HULL)  exit(-l) j 
if  ((inflle2=fopea(Bov«_iiame,"r">}==NULL)  axit(-l) ; 


for(tlLeta=angle_8t«Lrt;tlietu<=angl«_start-*-(orient_niiml  .0}*angle_tune'*-1.0; 
tlieta-*-=angle_tune)  { 
nvimber+sl ; 

sprintf  (buff  er_d,"X<i"  ,nui#ber) ; 
strcpy(buffer_a, Strings) ; 

Btrcpy(buffer_b, Strings) ; 
strcpy (buf f «r_c ,string4) ; 

8treat(buff er_a, buffer _d) ; 
strcat(buff er_b,buffer_d) ; 
streat (huff er_c ,buff er_d) ; 
strcat(buff er.a.stringl) ; 
streat (buf f er.b  ^strlngl) ; 
streat (buff «r_c .stringl) ; 


— -Gabor  Type  Kemal  Generator 
infile3=fopen("kem.dat",”«") ; 

ex=ODegaecoa(thetaepi/180} ; 
yy=ODega«sin(theta*pi/180) ; 
for  (i=index;l>=-index;i — ) 

( 

for  (j=index: j>=-index!j — ) 


} 


image=exp(-  (float)  (lei-t-J*  j)/(Bipaa*signa) ) : 

image l=cos(2.0epi* (vx* (float) i+vye (f loat ) j ) 4pba8e*pl/ 180 . 0) elmage ; 
lmage2=8in  (  2 . 0*pi*  (vx*  (float)  i+vy*  (float )  j )  'fphaseepi/ 180 . 0)  elmage ; 
rpartt [Index-1] [index- j]=lmagel ; 

rpartx[index-i]  [index- j]=-iBnge2*vi*2.0*pi-lBagel*(float)i*2.0/(8igmn*8igma) ; 
rparty [Index-i]  [index- j]=-inage2*vy*2.0*pi-imageie (float)  j*2.0/(8igBa*Bigma) ; 
fprintf (lnfile3,"%f\t",lfflagel) ; 


rprlntf (infile3,"\n") ; 

} 
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/* - convolution - */ 

If ((outl=fop«a(buffoK_a,''tf'*))»HULL)  oxitCO) ; 
il((ottt2=fop«n(btiff«_b,"w'*))MNULL)  axit(O) ; 
if  (  (out3=f op«n  (buf f or.c ,  "»*•)  )==iniLL)  oxlt  (0) ; 

for  (row=iad*x;row<k*m8ixo;row++)  /*  Initialixo  tb*  buffer  •/ 

{ 

for  (col=0;eol<CQLUMNSlZE;col++) 

•t 

fectoif  (infile  l.'*Xf"  .kbufferl  [row]  [col]) ; 
f  scenf  (inf ile2 .  "tt”  .*buf f er2  [row]  [col]  ) ; 
bufferl [row- index] [col] sbufferl [index] [col] ; 
buff er2 [row- index] [col] =buffer2 [index] [col] j 

> 

> 

for<i*0;i<R0WSIZE:i++)  /•  do  the  convolution  */ 

{ 

for  (  j=Oij<COLUMNSIZE;j++) 

imagelsO.O; 

iiB«ge2«0.0: 

iBege3=0.0; 

for  (rovl=-index;rowl<index+l;rowl++)  <  /*  convolution  unit  •/ 
for  <coll=-index:coll<lndex+l;coll++)  { 
col2=j+coll: 
if(col2<0}  col2s0; 

if (col2>COLl«NSIZE-l)  col2=CaLQMNSIZE-l ; 

image l+sbufferlbowl-findex]  [col2]*rpartt[iadox-*'rovl]  [index+eoll] ; 
iiiiage2-*-=buffer2[ravl-»indox]  [col2]«rpartx[indox'*’rovl]  [iadex-*-coll] ; 
iBago3-*'=buffer2[rovl'i’index]  [col2]  erparty  [indax+rovl]  [indox-*-coll] ; 

} 

> 

fprintf (out l,"Xf\t", image 1) ; 
fprintf  (out2/'Xf\t’',image2) ; 
fpr intf (out3 , "Xf \t " , imaged) ; 

)  /*  for  loop  for  J  */ 


fprintf (outl,"\n") ; 
fprintf (out2,"\n") ; 
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fprintf (outa,"\n") ; 

Hot  (eoloOicoKCOLUMNSIZEscol-M-)  f*  ahlftor  */ 

■C 

for  (rovsO ;  rov<k«mslsa  - 1 ;  rov44) 

< 

buff «rl  [rov]  [col]  =buf f orl  [rou+l]  [col]  ; 
buffers  [rov]  [col]  =buf f  orS  [row-*-!]  [col]  ; 

} 

if  (1<R0USIZE-Ind«x-1) 

{ 

f scanf  (inf  ilol .  "Xf  "  .tbuf f orl  Dcomsizo-l]  [col]  } ; 
fscanf  (liifll«3,"Xf"tkbuff«r3[k«raslz«-l]  [col]) ; 

> 

> 

> 

rovlnddnf  ilol) ; 
rovlnddnf  iloS) ; 
fclos«(outl}: 
f close (outs) ; 
f close (outa) : 

> 

fclosedafllel): 
f close (inf ileS); 

/* - Least  S<iuare  Estifflatlon - •/ 

if  (  orlent_num  «  S.O  )  solve3(R0WSIZE,C0LVNNSIZE.thzesold} ; 

else 

solve4(R0USIZE,C0LUI«SIZE.tliresold) ; 

> 


9  LEASTSQ.C 


/* 

*  Ibis  program  Is  to  minimize  tbe  motion  error  by  using 
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*  laast-aquara-tatlmatlon.  Tvo  fuaetloaa  art  laeludad: 

*  (1)  aolvadO :  function  for  rotating  4  orlantatlona  of  Gabor  typa  function 

*  (2)  8o1vo2():  function  for  rotating  2  orlantatlona  of  Gabor  typo  function. 

*  raqulrad:  ddd.dat,  tha  diffaranca  of  firat  fraaa  and  aacond  frama  of  iawgoa 
*/ 


tlncluda  <8tdlo.h> 
flncludo  <8trlng.b> 
dlncluda  <8tdlib.h> 
dlnclude  <math,h> 


aolva4(R0W, COLUMN,  throaold) 

Int  ROW. COLUMN; 
float  tbraaold; 

< 

int  l,j,k,lndax: 
float  yt[6],yx[6],yy[S]; 
float  u,v,ball; 

float  nuffll,num2,aua3,nuB4.num5: 
FILE  *ftl,*ft2,*ft3,*ft4,*ftS: 
FILE  afxl,afx2,afx3,afx4,afx5: 
FILE  '*fyl,*ly2,afy3,*fy4,*fy6; 
FILE  *outl,*out2,*out3,*out4; 


if  ((fxl=fopanC"xl.dat","r''))=*NULL>  axlt(O) ; 
if  ((fx2=fopan("x2.dat",'’r"))=NULL)  axltCO) ; 
if ((fx3"fopan("x3.dat","r"))®=HULL)  oxlt(O) ; 
iJ((fx4=lopan("x4.dat","r"»=NULL)  axitCO) ; 

if ((fyl=fopan("yl.dat"."r"))=«NULL)  axit(O) ; 
if ((fy2*fopan("y2.dat","r"))=®NULL)  axlt(O) ; 
if  ((fy3=fopan("y3.dat","r"))“NULL)  axit(O) ; 
if ((fy4=fopan(”y4.dat","r"))«»NULL)  axlt(O) ; 


if  ((ftl=fopaa<"tl.dat",''r"))=*NUU)  axit(O) ; 
if ((ft2=fopaa("t2.dat","r"))=aNULL)  axit(O) ; 
if  ((ft3=fopan("t3.daf',"r"))=HULL)  axit(O) ; 
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if  ((ft4=fop«ii{"t4.diif  ,"r"))— NULL)  •xit(O) ; 

If  ((out3sfop*n("dd<l.iag'',*'r‘’}}=BHULL}  MltCO) ; 
if ((outl=fop«n("u.dat","w"))«=HULL)  •xit(O) ; 

if  ((out2=fop*n(''v.dat",*'w"))=HULL)  axitCO) ; 

iad«x=0: 

for(i=0;i<R0W:i++) 

{ 

f  orC  jaO ;  j  <COHmN  I  j++) 

fscanf  (fxl,"Xf''  lAyxCO]) : 
f  scMif  (fx2 .  "Xf "  .*yx  Cl]  ) ; 
f  acanf  (f  x3 ,  "Xf  "  .*yx  [2]  > ; 
fscanf  (f  x4 ,  "Xf "  ,kyx  [3] ) ; 

f8canf(fyl,"Xf",*yy[0]); 
f8canf(fy2,"Xf",*yyCl]); 
fscanf(fy3."X«",tyyC2]); 
fscanf  (f  y4 ,  "Xf "  ,*yy  [3] ) ; 

fscanf (ft 1 . "Xf " ,*yt [0] ) ; 
fscanf (f t2 . "Xf " ,*yt [1] ) ; 
fscanf  (f  t3 ,  "Xf  "  ,*yt  [2]  ) ; 
fscanf (ft4 , "Xf " .*yt [3] ) ; 

fscanf (cuts , "Xf " ,*ball) ; 

numl=0.0; 

niui2=0.0; 

num3=0.0: 

nu]n4s0 . 0 ; 
numS=O.Oi 

if(  ball  s=  0.0  )  {  /*  datazmina  if  notion  occurs  a/ 

fprintf (out 1 , "Xf \t" ,0.0); 
fprintf (out2 , "Xf \t" ,0.0); 

} 

alsa  { 

fox  (k>0;k<4;ka4)  /a  /lO.O  fox  raducing  out  of  ranga  a/ 

{ 

nunl +=yx [k] /SO . Oayx [k] /60 . 0 ; 


SOFTWARE  LISTS 


49 


nuoS-i'syy  [k]  /SO .  0*yy  [k]  /SO .  0 ; 
nuia34=yy  [k]  /SO .  0*yx  [k]  /SO .  0 ; 

nuiB4+=yt  [k]  /SO ,  0*yx  [k]  /SO .  0 : 
num6+=yt [k] /SO . 0*yy [k] /SO . 0 ; 


if  (f  aba  (nviml*nuiii2-num3*nuin3)<=  thxasold) 

{ 

putsC  data  overflow  "); 

u=0.0; 

v=0.0: 

index+=l ; 

printf("\n  row=Xd,  coluiim=Xd\n'*,i,j); 

> 

else 

•C 

/•  Coordinate  System  would  be  Consistent  with  HATUlB  Coordinate  System  */ 
•s-  (num2*nuD4-nuffl3*nufflS)  /  (atiml*num2-nua3*nuD3) ; 
us-  (numl«numS'-nuD3*nun4)  /  (nual*nnm2-nua3*nuiB3} ; 

} 

fprlntf (outl , "Xf \t" ,u) ; 
fprintf (out2 , "Xf\t " , v) ; 

}  /*  ball  ==  0  */ 

y 

fprintf (outl, "\n") ; 
fprintf (out2 , "\n") ; 

} 

prlntfC  \ntatol  dataflow  number  is  Xd\n" , index) i 
removeC'tl.dat") ; 
remove("t2.dat") ; 
remove("t3,dat"} ; 
remove ("t4.dat'') ; 
removeC'xl.dat") ; 
remova("x2.dat") ; 
remove("x3.dat"} ; 
remove("x4.dat") ; 
removeC'yl.dat") ; 
remove("y2.dat") ; 
remove("y3.dat") ; 
remove ("y4.dat") ; 
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solves (ROW .COLUMN .thresold) 
int  ROW, COLUMN; 
float  tliresold; 

■t 

int  i.j.k, index; 

1  lost  yt  [2]  ,yx  [2]  ,yy  [2] ; 

float  numl.num2,aom3,nuB4,nua6.u.v,ball; 

FILE  eftl,*ft2; 

FILE  efxl,efx2; 

FILE  *fyl,efy2; 

FILE  *outl,eottt2,*oat3; 

if  ((fxl«fopen("xl.dat'',"r'’))«aMULL)  axit(O) ; 
if  ((fx2*fopen("x2.dat«,"r"))«NULL)  exit(O) ; 

if  ((fyl»fopen("yl.dat'',"r"))=«NULL)  exit(O) ; 
if  ((fy2=fopen("y2.dat","r"»«NUlL)  exit{0) ; 

if  ((ftl=fopen("tl.dat'’,"r"))«*NOLL)  exit(0> ; 
if  ((ft2*fopenC’t2.dat","r"))*»NULL)  exit(O) ; 
if  ((out3sfop«n("ddd.i]Bg",'*r"))ssMULL)  exit(O) ; 

if((outl=fopen("n.daf',"ii'’))=*NOU)  exit(O) ; 
if  (  (ont2--f  openCv-dat" ,  "w")  )  »NULL)  exit  (0) ; 

index^O ; 

for(isO;i<ROM;i++) 

{ 

fcr(j=0;j<C0LUMN;j++) 

i 

f8canf(fxl."Xf",fcyx[0l); 
f  scanf  (f  x2 ,  "Xf  "  ,Ryx[l]  ) ; 

f scanf (f yl , "Xf " .Ryy [0] ) ; 
f  scanf  (fy2 ,  "Xf"  .fcyytl] ) ; 
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f scanf (f tl , "%f " ,tyt [0] ) ; 
f  scanf  (f  t2 ,  ''%f”  ,kyt  [1] ) ; 
f  scanf  (cuts,  ''Xf  ‘  ,kball) ; 

nual=0.0; 
niim2=0 . 0 ; 
niun3=0.0; 
nuai4=0 . 0 ; 
iium5=0.0; 

if(  ball  =  0.0  )  {  /*  detarmine  If  motion  occurs  */ 

fprintf  (outl,"‘/{f\t",0.0) ; 
fprintf (out2,"Xf\t",0.0) ; 

} 

else  { 

for  (k=0;k<2ik++)  /•  /lO.O  for  reducing  out  of  range  •/ 

{ 

numl+=yx  [k] /BO .  0*yx  [k] /SO .  0 ; 
nuiB2+=yy  [k]  /BO .  0*yy  [k]  /SO .  0 ; 
nun3+=yy [k] /BO . 0*yx [k] /SO . 0 ; 
nuBi4+=yt  [k]  /BO .  0*yx[k3/S0 .0 ; 
nuiBB+=yt  [k]  /BO .  0*yj  [k]  /SO .  0 ; 

} 

if  (f  abs  (numlenum2-nuiii3*nuffl3)  <=t]iresold) 

puts("  data  overflov  "); 

u=0 . 0 ; 

v=0.0; 

index+=l: 

printf ( "\n  row=Xd ,  column=Xd\n" , i , j ) ; 

} 

else 

i 

t*  Coordinate  System  Transform  for  U  and  V  Otherwise  U  and  V  change  each  other*/ 
v=-  (nuffl2*n\iffl4~num3*nvu[£)  /  (nuBl*nuin2~nuffl3*num3) ; 
u=-  (numl  *numS-num3*num4)  /  (nunl  *num2-nuffl3*nuia33 
} 

fprintf (outl , "Xf\t" ,u) ; 
fprintf (out2 . "Xf \t" ,v) ; 

}  /♦  ball  ==  0  */ 

} 
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f pr intf (out 1 , "\n" ) ; 
fprintf (ottt2,"\n"): 


} 

printfC  \ntatol  dataflow  number  is  XdVn", Index); 

removeC'tl.dat") ; 

remove("t2.dat") ; 

removeC'xl  .dat") ; 

remove ("x2.dat") ; 

removeC'yl.dat") ; 

remove ("y2.dat'') ; 


> 


10  BMll.B 

studio  { 

from  -30  226  120  //set  up  position  of  the  camera  here  and  next  line 

at  -10  10  10 

up  0  0  1 

angle  27.1 

res  120  100 

aspect  1.2 

antlallas  adaptive 

background  (001) 

ambient  .8  .8  .8 

) 

/«  light  source  definition*/ 

light  {  type  point  falloff  1  position  60  120  80  color  26  25  26  > 

/*  objects  deflnatlon*/ 

//  for  red  x-axls  // 

surface  {  dlff  (.3  .3  .3}  shine  20  ,5  .5  .5  } 

cone  {  apex  000  base  26  0  0  apex.radius  1  base.radlus  1  } 
cone  {  apex  25  0  0  base  27  0  0  apex.radius  1  base.radlus  5  } 
cone  {  apex  27  0  0  base  35  0  0  apex.radius  5  base.radius  0  } 
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sphere  {  center  35  E  0  rsdlua  1  > 

sphere  {  center  42  S  0  radius  1  > 

sphere  {  center  35  5  10  radius  1  } 

sphere  {  center  42  5  10  radius  1  } 

cone  {  apex  35  5  0  apex.radius  1 
base  42  5  10  base_radius  1  } 
cone  {  apex  35  5  10  apex.radius  1 
base  42  5  0  base.radius  1  } 

//for  green  y-axis  // 

surface  <  diff  (0  1  0)  shine  20  .5  .5  .5  } 

cone  {  apex  000  base  0  25  0  apex.radius  1  base.radius  1  } 
cone  •(  apex  0  25  0  base  0  27  0  apex.radius  1  base.radius  S  } 

cone  {  apex  0  27  0  base  0  35  0  apex.radius  5  base.radius  0  } 

sphere  {  center  5  36  10  radius  1  > 
sphere  <  center  12  35  10  radius  1  > 
sphere  {  center  8.5  35  0  radius  1  > 

cone  {  apex  8.5  35  5  base  5  35  10  apex.radius  1  base.radius  1  } 
cone  {  apex  8.5  35  5  base  12  35  10  apex.radius  1  base.radius  1  } 

cone  {  apex  8.5  35  5  base  8.5  35  0  apex.radius  1  base, radius  l  } 

//  for  blue  z-axls  // 

surface  {  diff  (0  0  1)  shine  20  .5  .5  .5  } 

cone  {  apex  000  base  0  0  25  apex.radius  1  base.radius  1  } 
cone  {  apex  0  0  25  base  0  0  27  apex.radius  1  base.radius  6  > 

cone  {  apex  0  0  27  base  0  0  35  apex.radius  6  base.radius  0  } 

sphere  {  center  -5  -5  25  radius  1  > 
sphere  {  center  -12  -5  25  radius  1  > 
sphere  {  center  -5  -5  35  radius  1  > 
sphere  {  center  -12  -5  35  radius  1  > 

cone  {  apex  -12  -5  35  base  -5  -5  25  apex.radius  1  base.radius  1  } 
cone  {  apex  -5  -5  35  base  -12  -5  35  apex.radius  1  base.radius  1  > 

cone  {  apex  -5  -6  26  base  -12  -6  25  apex.radius  1  base.radius  1  > 


/*  Sphere  */ 
surf  { 
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diff  (  .8  .1  .1) 
spec  .8  .1  .1 
shlse  30 
} 

sphere  {  center  10  60  10  radius  10>  //  set  up  position  of  sphere  here 

/*  The  end  •/ 


r 


